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PREFACE 


The  report  for  the  work  performed  under  Contract  No.  D.O.  0064, 
DAAH01— 87— D-0021  consists  of  three  volumes:  (1)  "Extension  and  Updating  of  the 
Computer  Simulation  of  Range  Relative  Doppler  Processing  and  Invariant  Mapping  for  MM 
Wave  Seekers",  (2)  "User  Manual  of  the  Range  Relative  Doppler  Processing  Simulation  for 
MM  Wave  Seekers",  and  (3)  "Computer  Simulation  of  a  Doppler  Beam  Sharpening  Radar  for 
MM  Wave  Seekers".  The  period  of  performance  is  April  12  to  December  31,  1988. 

The  main  objective  of  Volume  I  is  to  extend  and  update  the  MM  wave  computer 
simulation  developed  for  Contract  No.  DAAH01— 87— D— 0021,  D.O.  18,  documented  in  the 
final  report  dated  February  1988,  and  entitled  "Increasing  Azimuth  and  Elevation  Resolution  of 
MM  Wave  Seeker  Systems  Using  Coherent  or  Noncoherent  Range  Relative  Doppler 
Processing  (RRDP)  with  Constant  or  Linear  Frequency  Modulation  and  Invariant  Mapping". 

Volume  II  is  a  User  Manual  for  the  computer  simulation  documented  in  Volume  I. 
With  this  manual,  MICOM  personnel  should  be  able  to:  (1)  install  the  Fortran  software  on  an 
IBM  PC/compatible,  (2)  duplicate  the  results  presented  in  this  report,  and  (3)  run  the 
simulation  for  different  clutter  maps  and  targets. 

The  main  objective  of  Volume  III  is  to  develop  Doppler  Beam  Sharpening  for  MM 
wave  seekers.  In  this  application  the  geometry  is  significantly  different  That  is,  the  seeker 
follows  a  straight  horizontal  trajectory  while  the  antenna  performs  a  forward  near  r  ricular  scan 
perpendicular  to  the  trajectory. 


Robert  J.  Polge 


Huntsville,  Alabama 
January  1989 
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L  Introduction  and  Background 

One  of  the  most  serious  deficiencies  in  the  existing  MM  wave  seeker  systems  is  the 
lack  of  azimuth  resolution.  As  a  result,  detection  and  discrimination  of  targets  in  a  clutter 
background  is  quite  difficult.  Azimuth  resolution  cannot  be  increased  using  standard  SAR 
techniques  because  the  MM  wave  nonlinear  geometry  magnifies  the  problems  of  range  walk 
and  doppler  bias. 

Two  new  digital  processing  techniques  were  developed  [1,2]  to  extend  the  application 
of  SAR  processing  to  nonlinear  trajectories.  They  are  called  "Range  Relative  Doppler 
Processing"  and  "Invariant  Mapping".  SAR  Doppler  Processing  of  each  range  cell  was 
extended  to  nonlinear  trajectories,  by  using  as  Doppler  reference  a  fictitious  scatterer  at  the 
center  of  the  range  cell.  This  resolves  the  problems  of  range  walk  and  doppler  bias.  In  [1]  and 
part  I  of  [2]  the  approach  was  to  perform  a  time  Taylor  expansion  with  respect  to  the  fictitious 
scatterer,  the  effect  of  azimuth  on  doppler  being  introduced  later  as  a  parabolic  fit  at  the  center 
and  extremities  of  the  range  cell.  The  simulated  return  of  the  fictitious  scatterer  is  used  for 
heterodyning.  After  deramping,  windowing,  and  FFT  processing,  the  range  cell  is  divided  into 
trapezoidal  azimuth  cells  corresponding  to  a  constant  frequency  increment  Af,  where  Af  is  the 
inverse  of  the  observation  interval,  DQjJ. 

Pan  II  of  [2]  extended  the  application  of  Range  Relative  Doppler  Processing  (RRDP)  to 

linear  frequency  modulation  (discrete  or  continuous).  Another  important  contribution  was  that 

the  round  trip  delay  model,  and  therefore  the  signal  model,  includes  explicitly  not  only  time 

but  also  azimuth  and  elevation  increments.  The  azimuth  and  elevation  angular  increments 

(|ij,e-)  for  an  arbitrary  scatterer  in  the  jth  range  cell  are  measured  with  respect  to  the 

■** 

reference  fictitious  scatterer  C j .  It  is  clear  that  both  the  signal  return  and  the  doppler  due  to 
-* 

C-  are  functions  of  (|-L.,e-),  where  the  effect  of  £•  is  small  in  the  case  of  a  narrow  range  bin. 
The  round  trip  delay  r.  to  C  is  a  function  of  (t,p-,£|)  where  t  is  relative  time.  A  good 
approximation  for  the  round  trip  delay  Tj(t,Pj,Ej)  is  obtained  using  a  three-dimensional  Taylor 
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expansion  about  the  fictitious  scatterer  state.  As  before,  the  fictitious  scatterer  at  the  center  of 
the  range  cell  serves  as  zero  doppler  reference.  This  means  that  the  return  signal  must  be 
heterodyned  with  a  synthetic  signal  simulating  the  return  of  the  fictitious  scatterer,  including 
the  effect  of  frequency  modulation  when  applicable.  The  remainder  of  the  processing  is 
unchanged:  (1)  A/D  conversion,  (2)  deramping,  (3)  windowing,  and  (4)  FFT  processing.  It 
was  also  shown  that  the  stability  requirement  during  the  observation  interval  can  be  expressed 
as  a  bound  on  frequency  drift  which  is  practically  attainable. 

This  report  extends  and  updates  the  Range  Relative  Doppler  Processing  (RRDP) 
presented  in  [2],  and  also  replaces  the  invariant  mapping  presented  in  [1]  by  a  new  and  a  much 
more  efficient  algorithm.  Frequency  modulation  is  not  included  in  the  simulation  at  this  time. 
Most  of  our  efforts  on  the  RRDP  simulation  were  concentrated  on:  (1)  a  more  general 
geometry  where  the  trajectory  is  not  restricted  to  a  vertical  plane,  (2)  extending  the  RRDP 
simulation  to  include  the  effect  of  scatterer  height,  (3)  entering  interactively  the  parameters 
and  the  clutter  data  which  means  that  only  the  changes  need  to  be  typed  online,  (4)  an  exact 
simulation  for  the  dish  antenna  instead  of  a  sin  x  /x  pattern,  (5)  a  reduction  of  impulse 
invariance  over  the  entire  footprint  through  frequency  interpolation  and  a  nearly  perfect 
compensation  for  range  and  antenna  gain.  All  the  subroutines  in  the  RRDP  simulation  were 
thoroughly  tested  and  most  of  them  were  significantly  improved. 

Much  efforts  were  expanded  to  develop  a  better  and  much  more  efficient  invariant 
mapping  technique.  The  heart  of  this  technique  is  an  algorithm  which  computes  the  area  of 
c^erlap  between  a  rectangle  of  the  absolute  grid  and  a  trapezoidal  azimuth  cell  from  the 
current  range  cell.  The  mapping  is  area  invariant  which  means  that  if  a  trapezoid  overlaps 
with  three  rectangles  the  sum  of  the  three  overlaps  is  equal  to  the  area  of  the  trapezoid. 

Using  SURFER  from  Golden  Software  Inc.  [7],  we  can  display  three-dimensional 
(3-D)  maps  or  contour  maps  before  and  after  mapping  on  our  Compaq  386.  In  the  previous 
report  we  were  using  Univac  DISSPLA  which  is  a  much  slower  procedure.  For  two— 
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dimensional  plots  we  use  GRAPHER  which  is  also  from  Golden  Software. 

The  organization  is  as  follows.  The  geometry  and  dish  antenna  simulation  are 
discussed.  Then,  an  extended  Taylor  expansion  for  the  round  trip  delay  x.(h,t,p.,e),  which 
includes  scatterer  height,  is  presented.  Two  simulation  programs  are  developed.  The  first 
program,  FRRDP  generates  a  data  file  FRRDPG.GRD  which  stores  reflectivity  versus  azimuth 
and  range.  The  second  program  CRRDP  uses  the  new  invariant  technique  to  generate  data 
file  CRRDPG.GRD  which  stores  reflectivity  versus  absolute  x— y.  The  graphics  software 
package  SURFER  transforms  the  reflectivity  data  files,  FRRDPG.GRD  and  CRRDPG.GRD, 
into  3— D  displays  or  contour  plots.  The  simulation  is  illustrated  by  examples.  This  report  is 
concerned  with  the  theory,  implementation  and  performance  of  the  simulation.  An  User 
Manual  is  also  available  [8]. 

2.  Geometry  and  System  Parameters 

Figure  2.1  shows  the  geometry  for  a  footprint  at  central  time  t  ,  where  q(t  )  defines  the 
center  of  the  footprint.  In  our  previous  work  the  seeker  was  assumed  to  remain  in  a  vertical 
plane,  this  limitation  is  now  removed.  The  only  assumption  left  is  that  the  vector  velocity 
remains  constant  within  the  observation  interval  D  ^ 

'V^VV1!1  '  <21) 

-* 

The  vector  m  defines  the  position  of  the  antenna  at  central  time,  and  its  rectangular 
coordinates  with  respect  to  the  absolute  origin  6  are 

m  =  a(tc)  =  {x(m)  =  dx,  y(m)  =  dy>  h).  (2.2) 

Figure  2.1  shows  an  absolute  system  of  reference  {0x,Oy,Oz}  where  6  is  a  fixed  point 
on  the  ground.  The  positions  of  the  seeker,  the  current  footprint  center,  and  the  scatterers  are 
normally  defined  in  absolute  coordinates.  For  example,  the  footprint  center  is  given  by 

5*  =  q(tc)  =  {xq,yq,°}.  (2.3) 
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However,  within  one  observation  interval,  a  spherical  representation  of  a  ground  point  with 
respect  to  the  central  seeker  position  m  is  more  useful:  {range,  azimuth,  elevation},  where 
azimuth  is  referred  to  Oy  and  elevation  is  referred  to  the  downward  vertical.  For  example,  the 
footprint  center  in  the  spherical  representation  is  given  by 

q*-  m  =  {rq,  aq=  a*=  Angle(q*-mg,  Oy),  (3q=  (5  =  Angle(mg-m,q  -m)}.  (2.4) 

Similarly,  let  C-  be  the  ith  scatterer  and  let  be  its  ground  projection,  then 

Cgi-m=  {fgi.ctj.  Pi)  .  (2.5) 

In  the  remainder  of  this  report  the  azimuth  and  elevation  angles  of  a  ground  point  are  always 
with  reference  to  m. 

The  positions  of  the  antenna  within  Dq^,  and  of  a  scatterer  C.  located  hr  above  the 
ground  are  easily  written  in  absolute  coordinates. 


a(ta)  =  a(t  +  tc)  =  {dx  +  vxt,  dy  +  vyt,  h  +  vzt} 
D  u  D, 

-°S  <;!<_£ 


(2.6a) 


(2.6b) 


where  t  denotes  absolute  time,  and  t  is  the  relative  time.  The  absolute  coordinates  of  the  ith 

cl 

scatterer  C  are 

Cj  =  {xp  yj,  fy}  .  (2.7) 

It  will  be  shown  that  the  height  of  the  scatterer  h|  significantly  affects  detection.  The  ground 
projection  of  the  scatterer  is  defined  in  absolute  coordinates  as 

egi  =  txi*  yi- 01  ■  <2-8) 

It  is  obvious  that  a  scatterer  above  the  ground  can  be  read  as  the  sum  of  its  ground  projection 
and  a  vertical  vector 
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ci  =  cgi  +  (0, 0,  \) 


(2.9a) 


Ci  =  (C  .  -  m)  +  (0,  0,  fij)  +  m  (2.9b) 

where  (2.5)  or  (2.8)  are  used  depending  on  the  preferred  system  of  reference.  Assuming  that 

•4  4 

the  scatterer  projection  C  •  is  somewhere  within  the  jth  range  cell,  then  C-  may  actually 

appear  in  the  (j  -  l)th  range  cell,  if  the  height  hj  is  significant. 

Within  one  range  cell  the  reference  is  a  fictitious  scatterer  at  the  center  of  the  range 

-» * 

cell.  Figure  2.2  shows  the  jth  range  cell  where  the  reference  Cj  is  the  flat  fictitious  scatterer 

•4 

at  the  center  of  the  cell.  This  scatterer  is  defined  with  respect  to  m  by 

C*  —  m  =  {r*,  a*,  ($* }  ,  (2.10) 


*  *  -•*  -» * 
where  a-  =  a  because  C-  is  on  the  ground  line  of  sight.  Since  we  want  to  use  C-  as  a 
J  J  J 

reference  for  the  jth  range  cell,  we  define  the  azimuth  and  elevation  of  C  •  in  terms  of  small 

O1 

* 


incremental  variations  (p^  ep  about  (<Xj ,  ) 


ai  =  a  +Hi 


(2.11a) 

(2.11b) 


Give.i  the  position  of  the  jth  fictitious  scatterer,  the  position  of  the  ground  projection  of  an 


arbitrary  scatterer  within  the  jth  range  cell  is  entirely  defined  by  the  angles  p-  and  e-,  which 

will  be  denoted  as  relative  azimuth  and  relative  elevation.  The  position  of  C-  follows  from 

that  of  the  ground  point  C  ■  by  an  upward  vertical  motion  h,. 

As  before,  the  footprint  is  defined  as  the  3dB  contour  illuminated  by  the  antenna  at 

central  time.  While  the  footprint  moves  within  the  observation  interval,  most  of  this  motion  is 

compensated  by  Range  Relative  Doppler  Processing  (RRDP).  Given  the  antenna  position  m 

and  using  (2.3)  or  (2.4),  the  footprint  is  entirely  defined  by  its  center,  that  is  {x  ,y  )  or 
*  *  .....  99 

{a  ,p  }.  The  3dB  footprint  is  divided  into  range  cells  of  width  Ar,  and  for  each  range  cell, 
the  median  range  and  the  interval  for  p  {-pmx  <  p  <  pmx}  must  be  computed.  Within  the  jth 
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z 


Figure  2.2  Fictitious  scatterer  and  ith  scatterer  in  the 
jth  range  cell: a.  -  a  +  Fj . 


Y 
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* 

range  cell,  c  is  zero  at  medium  range  and  maximum  at  minimum  range  (r  •  —  Ar/2).  The 

J 

maximum  values  of  (e,n),  that  is  (emx,  /xmx)  are  used  to  determine  which  terms  in  the 
Taylor  series  expansion  of  round  trip  delay  can  be  neglected.  Appendix  A  contains  all  the 
formulas  which  define  the  relevant  information  in  a  footprint. 

The  ensemble  of  scatterers  is  denoted  as  clutter  map.  Details  of  the  clutter  map 
simulation  are  in  section  7.  There  are  two  options  to  enter  the  ith  scatterer  in  the 
simulation:  {x.,y.,h-}  or  For  each  scatterer  C-,  defined  in  rectangular  or 

spherical  coordinates,  one  needs  to  compute  the  following:  (1)  r.,  range  to  the  scatterer, 
(2)  the  index  of  the  range  cell  containing  C^,  (3)  the  elevation,  azimuth  and  their  relative 
counterparts  (a;,  /3-,  /x-,  e;),  and  (4)  the  antenna  gain  in  the  direction  (C.-4(t  ))  at  the 
start,  middle,  and  end  of  the  observation  interval.  The  antenna  gain  at  an  arbitrary  time 
within  the  observation  interval  follows  by  parabolic  interpolation.  This  information 
together  with  the  Taylor  coefficients  computed  in  Appendix  D  is  all  that  is  needed  to 
simulate  the  signal.  Appendix  B  contains  all  the  formulas  which  define  the  relevant 
information  for  a  scatterer  C-. 

In  this  report  the  seeker  transmits  simple  RF  pulses  without  linear  frequency 
modulation.  The  system  parameters  are  given  in  Table  2.1.  The  system  operates  in  the 
forward  looking  SAR  strip  mode. 

3.  Dish  Antenna  Simulation 

In  our  previous  work  the  dish  gain  pattern  was  approximated  as  a  sin  x  /x  function. 
An  exact  dish  gain  pattern  is  now  implemented.  The  gain  pattern  for  a  circular  dish  with 
diameter  D  at  direction  sine  s  =  sin/3  is  given  by  [3,4], 

J,  (  7TS  D/ A) 

G(sD/A)  =  2A  -j-wsY)/\)  (31) 

where  is  the  Bessel  function  of  the  first  kind  of  the  first  order,  and  A  is  the  area  of  the 
circular  dish.  Figure  3.1  shows  the  normalized  two-way  dish  pattern  Gn(w)  in  dB  where  w 
=  512sD/A. 
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Table  2.1:  List  of  System  Parameters 


Parameter 

Symbol 

Value 

Antenna  Beamwidth 

e 

43  mrad 

Frequency 

f 

94  GHz 

Range  Bin 

Ar 

1  m 

Total  Transmitter  Bandwidth 

B 

500  MHz 

PRF  or  FM  sweep  Repetition  Rate 

50  KHz 

CFAR  13  Cell  Parametric 

Cell-averaging  Line  CFAR 

Not  implemented  yet 

Log  Receiver 

Assumed 

ideal 

Seeker  Velocity  (input) 

a(tc)=qfc={vx,vy,vz)w:r.t.  0 
=  {av,(3v,v}w.r.L  in 

Seeker  Depression  angle 

Sv 

20° 

or  Slant  Angle  (inputs) 

Pv 

70° 

Seeker  Azimuth  (input) 

av 

35° 

Dwell  time  About  Central 

n  , 

20  ms 

Footprint  (input) 

Mob 

Position  of  Seeker  at  Central 
time  (input) 

a(tc)=m= 

{dx,dy,h } 

h  =  900  m 

Center  of  Footprint  (input) 

q(tc)=^* = { yq,0 }  w.r. t.  0 
={rq,a*,|3*}w.r.t.  m 

9 


10 


BEE 


The  one-way  dish  gain  is  computed  using  the  function  DISH(SINB,  EDDISH,  DDISH) 
where  SINB  =  sinP  =  the  off  boresight  angle,  IDDISH  =  { 1,2,3}  =  (rectangular,  circular  Dish, 
sin  x  /x},  and  DDISH  =  D fk  =  dish  diameter  in  X  units.  For  a  be  am  width  9  =  43  mrad,  the 
two-way  3dB  footprint  corresponds  to  9^  =  21.5  mrad.  One  can  compute  the  equivalent 
diameter  for  IDDISH  =  2  and  3  by  solving 

40  log10[DISH(6H,  2,  D2)]  =  -3  (3.3a) 

40  log10[DISH(6H,  3,  D3)]  =  -3.  (3.3b) 

The  results  are:  D2  =  17.15094  and  D^  =  14.80884.  Note  that  the  one-way  power  gain  is 
proportional  to  the  modulus  square  of  the  electric  field.  The  listings  for  the  Function  DISH 
and  the  associated  Bessel  function  BESSJ1  are  in  Appendix  C. 

4.  Taylor  Expansion  of  Round  Trip  Delay  T- 

Assume  that  the  ground  projection  of  the  ith  scatterer  within  the  jth  cell  is  specified  by 

-  m  =  (rgi,  a.  =  a*  +  Pj,  ^  =  (3*  +  .  (4.1) 

-4 

Appendix  B  shows  how  to  compute  all  the  relevant  information  for  the  scatterer  C,  and  in 

■4  4 

particular  r-,  e-,  and  p-.  From  Fig.  2.1,  the  elevation  of  with  respect  to  m  is 

Pi  =  cos  ^h/tgj)  •  (4.2) 

It  follows  that  the  absolute  coordinates  are  computed  by  combining  (2.9b)  and  (4.1). 

Table  4.1  gives  the  rectangular  coordinates  of  the  vectors  2'(t_)  =  £  '(t  ),  2(t  ), 

^  d  C  d 

■4  -4 

C,  G  -  £t(t).  One  can  now  compute  r,  the  slant  range  between  the  seeker  at  time  t  =  t  +  t 

11a  1  a  C 

-♦ 

and  the  scatterer  C, 

r;  |Cj  -  a(ta)|2  (4.3) 

After  few  manipulations  one  obtains 

r{  =  |Cj  -  3(ta)  j  =  [h2tan2Pi  +  v2t2  +  H2 

-2  h  tanp.  t  (v  sina.  +  v  cosa-}  +  2H  v  t]^2,  (4.4) 

i  x  i  y  i  z 


ii 


where  H  =  h  -  a.  =  a  +  ji^,  +  £.. 

During  RRDP  of  the  jth  range  cell,  the  reference  state  is  defined  by 

{m  =  a(tc),q*,C*}  ,  (4.5) 

and  the  reference  slant  range  is 

r*  =  |  C*  -  m  |  .  (4.6) 

The  slant  range  r-  can  be  expressed  in  terms  of  incremental  variations  about  the  reference 
range 

ri  =  r-flkj,  t,  ^ ,  ei  ;  m,  q*,  C*)  (4.7a) 

r*  =  r(0,0,0,0  ;  m,  q*  C*)  .  (4.7b) 


Table  4.1.  Rectangular  coordinates  for  relevant  vectors 


Vector 

x— coordinate 

y-coordinate 

z— coordinate 

rty 

vx=vsinpy  sinay  vy=vsinPy  cosay 

vz=  — vcos(5y 

4(0 

d„  +  v  t 

d  +  v  t 

h  +  v  t 

a 

X  X 

y  y 

z 

ci 

h  tanfy  sina-+dx  h  tan^cosct+dy 

\ 

h  tanP-sina.-v  t  h  tan(3.cosa.-v  t 

i  1a  1  I  y 

— (h— hj+vzt) 

The  slant  range  r.  can  be  expressed  as  a  Taylor  expansion  in  {fr,,  t,  p.,  e-)  about  the 

reference  slant  range  r. .  Equivalently  one  can  expand  the  round  trip  delay  x-  =  =—», 

J  1C 

where  c  is  the  speed  of  light 

Since  the  formula  for  r.  is  complicated,  we  need  to  define  a  simple  general 
notation.  Assume  that  x=x^  is  a  function  of  four  variables  (v^,  yy  v4)-  The  partial 
derivative  of  x  with  respect  to  v^  evaluated  at  the  reference  state  is  denoted  as  x^.  This 
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notation  is  extended  to  multiple  partial  derivatives.  For  example 


x123"  dv1dv^Jv^  X 


vl=v2  =  v3  =  v4  =  ° 


(4.8) 


C-m| 


1  1  —  1 

Using  this  notation,  the  round  trip  delay  x;  =  2  — - —  becomes 

1  V 

*2  2 

x(  1,2,3, 4)  =  t+(x1Tt+t2t+x4e)  +  {tnr  +  x{2ht  +  x14te  +  x22  2  +  *23^ 


2  ^ 

+  X24te  +  X44  2  }+{xlll5  +x112  V+ X114  V  +  'C122^2_  +  t123ht^+X124r,ie 


h2e 


2_  3_  2_  2_  2  _  2 

+  x144  “2 f  +  X222  r  +  x223  2  +  X224  ¥  +  X233  +  X234  tJI£  +  X244  ^2" 


+x 


444  5" 


(4.9) 


where  the  subscript  i  has  been  dropped  to  further  simplify  the  notation.  The  coefficients 
of  the  Taylor  series  are  computed  in  Appendix  D.  The  Fortran  implementation. 
Subroutine  TAYLOR,  is  in  Appendix  E. 

Numerical  evaluation  of  the  coefficients  of  (4.9)  shows  that  many  of  them  are 
completely  negligible.  Therefore,  the  final  practical  formula  is 

x(h,t,n,e)  =  (x  +  Xjfi  +  x4e)  +  t(x2  +  x23n  +  x24e 
2  2 

+  *12^  +  *233  *¥  +  1  (i22  +  *223^  +  *224e)'  (410) 

The  accuracy  of  (4.10)  was  validated  by  comparing  the  simulated  round  trip  delay  and 
doppler  to  the  exact  values.  For  example,  the  simulated  doppler  is  (-6994.76  rad/sec) 
while  the  actual  value  is  (—6992.87  rad/sec).  It  follows  that  (4.10)  can  be  used  for  an 
accurate  signal  synthesis. 
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S.  Range  Relative  Doppler  Processing 

The  concept  of  Range  Relative  Doppler  Processing  (RRDP)  is  presented  in  Part 
I  of  [2]  and  the  implementation  for  linear  frequency  modulation  is  part  n  of  [2].  For 
completeness  a  concise  derivation  of  RRDP  follows  for  the  RF  case. 

The  return  corresponding  to  the  scatterer  C.  located  within  the  jth  range  cell  is 
Sj(t)  =  Bjft)  cos  (coQ(t  -  xp  +  §)  ,  (5.1) 

where  is  the  transmitter  phase,  x.  the  round  trip  delay  is  given  in  (4.10),  and  B.(t) 
incorporates  antenna  gain,  range  attenuation  and  scatterer  strength.  Using  the  Taylor 
expansion  (4. 10)  one  can  write 


£»(t  -  xp  +  %  =  coo(l  -  x2)t  +  y* 

Vi  =  -<oo[t(x23H  +  x24e  +  x12h  +  x233  ^) 

t2  _  .  _  _ 

+  1  {x22  +  W  +  x224e)  +  <x  +  xlh +  x4e)]  +  *• 


(5.2a) 


(5.2b) 


As  illustrated  in  Fig.  5.1  the  first  step  of  the  processing  is  mixing  at  the 
frequency  of  the  fictitious  scatterer. 


%j  *  %(>  -  V  •  (5-3) 

where  co  •  is  a  function  (nearly  linear)  of  j  through  x7.  The  synthetic  signal  exp(jco  t) 

yj  Z  OJ 

* 

is  the  estimated  return  from  Cj ,  a  fictitious  scatterer  at  the  center  of  the  jth  range  cell. 
This  means  that  the  analog  quadrature  components  are  given  by 

xdi  =  Bi(t)  cos  (5  4a) 


xqi  =  Bi(t)  sin  (5-4b) 

After  filtering  and  A/D  conversion,  deramping  is  performed  on  the  quadrature 

2 

components  as  a  matrix  rotation  to  cancel  the  quadratic  phase  a)Q  t  12.  The  deramped 
quadrature  components  are 
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cos(<0n.t) 


L5 


xdi  =  Bi(t)  cos 

(5.5a) 

xqi  =  Bi(t)  sin  (y(t)) 

(5.5b) 

where 

* 

y(t)  =  Vi  +  %X22  r 

2 

=  +  *24e  +  i  I2h  *  *233  2> 

,2  .  .  _  . 

+  ^  (^223^  +  x224e)  +  (x  +  t  j  h  +  x4e)]  +  %  (5.5c) 

and  the  residual  quadratic  phase  term  is  quite  small. 

It  follows  from  (5.5c)  that  the  relative  doppler  frequency  corresponding  to  the 

C-  scatterer  is  (neglecting  a  small  linear  frequency  variation) 

-  2 
<*>di  =  -  -“0^23**  +  *24£  +  *12h  +  *233  V1'  (5'6) 

When  (ft  =  0,  e  =  0),  i.e.,  on  the  median  of  the  range  cell,  (5.6)  becomes 

co^h  =  0,  e  =  0)  =  ~ <»0[^23m'  +  ^233  (5J) 

2 

where  the  |i  term  is  relatively  small. 

Figure  5.1  shows  that  the  processing  of  the  array  of  deramped  quadrature 
components  consists  of  three  steps:  windowing,  FFT,  and  scaling.  The  Kaiser  window 
is  recommended  [5,6]  because  it  allows  a  continuous  tradeoff  between  resolution  and 
leakage.  The  output  of  the  FFT  should  be  in  the  modulus/phase  format.  The  modulus 
must  be  scaled  to  compensate  for  antenna  pattern  and  range  attenuation.  The  FFT  cells 
are  numbered  (k=0,  N-l)  and  the  kth  cell  corresponds  to  a  relative  doppler  interval: 

(k-0.5)Ao)  <  o)^  <  (k+0.5)Aco  (5.8a) 

Aco  =  2k  Af  =  27t/DQb  (5.8b) 

where  DQb  is  the  duration  of  the  observation  interval  in  seconds. 
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A  point  scatterer  will  affect  three  adjacent  FFT  cells  as  illustrated  in  Fig.  5.2  for 
a  rectangular  window  and  for  a  Kaiser  window  with  parameter  equal  to  n.  Assume  for 
example  that  the  relative  doppler  for  a  point  scatterer  is  kAco,  this  means  that  the 
position  of  the  scatterer  in  the  jth  range  cell  is  defined  by  the  angular  increments  (e,(i) 
such  that  =  kAco  in  (5.6).  As  a  first  approximation  one  could  use  the  scaled 
modulus  A^  of  the  kth  cell  (on— target)  as  an  estimate  of  the  strength  of  the  scatterer. 
However,  A^  varies  with  the  relative  location  of  the  scatterer  within  the  range  cell, 
being  maximum  when  the  scatterer  is  located  at  the  center.  The  normalized  values  of 
the  (early,  on— target,  late)  components  are  shown  in  Fig.  5.3  for  two  cases:  (a) 
scatterer  at  the  center  of  the  cell,  and  (b)  scatterer  on  the  border  between  on— target  and 
late.  Theoretically,  one  could  use  an  interpolation  algorithm  to  define  the  center  of  the 
scatterer  more  precisely  from  the  (early,  on— target,  late)  components;  very  much  as  for 
range  interpolation.  However,  this  is  not  practical  in  the  case  of  multiple  scatterers. 
Instead  we  will  detect  a  single  point  scatterer  in  a  range  cell  as  three  scatterers  located 
at  the  centers  of  adjacent  cells.  Continuing  our  example  the  three  scatterers  will  be 
detected  at 

{(e=0,  nk—1)  ,  (£=0,  nk)  ,  (e=0,  nk+1)}  (5.9) 

where  (i^  is  the  solution  of  (5.7)  for  {co^  =  f  Ao),  t  =  k-1,  k+1).  The  value  of  p 
which  defines  the  border  line  between  the  kth  and  the  (k+l)th  FFT  cells  is  obtained  by 
solving  (5.7)  for  co^  =  (k+0.5)Ad). 

The  maximum  values  of  (£,p)  are  denoted  as  emx  and  pmx,  respectively. 
Appendix  D  explains  the  computation  of  (Emx,  pmx). 

Since  a  point  scatterer  will  also  affect  the  adjacent  cells  it  will  be  detected  as 
nine  samples  altogether.  In  a  typical  application  Emx  is  much  smaller  than  pmx  and 
formula  (5.7)  can  be  used.  It  indicates  a  nearly  linear  correspondence  between  doppler 
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rigure  5.2.  Normalized  response  for  both  rectangular  and  Kaiser  windows, 
Kaiser  (early, on, late)  samples  are  shown.  Central  target^,  □  a). 

Border-line  target  kg.*). 


42.9515 


p=0 


32.4764 


=  64 


Figure  5.3  FFT  line  spectrum  for  a  unit  sinewave  and  Kaiser  window((3=7i). 
N=64,f0=(m+p)Af  (p=0:central  target,p=0.5:borderline  target). 


2  - 

boundaries  and  angular  boundaries  because  (t  u.  /2)  is  much  smaller  than  (t  (i)  for 

kM’ 

the  useful  range  of  \i.  Yet  including  the  quadratic  term  increases  the  accuracy  by  about 
five  percent. 

6.  Reduction  of  Impulse  Invariance  Over  a  Footprint 

A  unit  scatterer  is  detected  as  a  narrow  three-dimensional  Gaussian  shaped 

pulse  [1,2].  The  impulse  corresponding  to  a  unit  scatterer  located  at  q  (the  center  of 

the  footprint)  will  be  defined  as  the  impulse  response  of  the  system.  Impulse 

invariance  means  that  a  unit  scatterer  anywhere  within  the  footprint  will  produce  the 

same  impulse  response.  It  is  difficult  to  achieve  impulse  invariance  over  the  entire 

footprint  because  of  range  attenuation,  antenna  gain  variation,  and  FFT  sampling 

effects.  The  results  in  [2]  show  the  impulse  response  varying  by  as  much  as  31%,  and  it 

was  recommended  to  correct  this  problem.  One  of  the  major  accomplishments  of  the 

present  effort  is  to  reduce  impulse  invariance  over  an  entire  footprint  to  within  ±4%. 

The  causes  of  impulse  variation  and  the  proposed  solutions  are  discussed  below. 

The  range  to  a  scatterer  varies  very  little  during  the  observation  interval,  but  it 

depends  significantly  on  the  location  of  the  scatterer  within  the  footprint.  The  range 

variation  is  even  greater  from  one  footprint  to  the  other.  Since  the  radar  equation 

involves  the  fourth  power  of  range,  the  range  attenuation  must  be  compensated  with 

respect  to  a  reference  range,  on  the  ground  map.  More  precisely,  the  range 

compensation  of  the  jth  cell  is  performed  by  multiplying  the  amplitude  of  the  FFT  cells 
*  4  * 

by  (rj/rref)  i  where  Tj  is  the  median  range.  This  compensation  works  very  well 
because  range  is  precisely  known  at  central  time  and  changes  very  little  during  the 
observation  interval. 

The  antenna  gain  varies  by  as  much  as  3dB  within  the  footprint.  It  follows  that 
a  scatterer  on  the  edge  of  the  footprint  will  produce  an  impulse  response  half  as  small 
as  a  scatterer  at  the  center  of  the  footprint.  For  each  FFT  cell  the  antenna  gain  can  be 
computed  at  central  time.  Therefore,  most  of  the  variation  due  to  antenna  gain  can  be 
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removed  by  dividing  each  cell  return  by  the  estimated  antenna  gain  for  that  FFT  cell. 

We  have  explained  that  the  detection  of  scatterers  was  transformed  into  a 
problem  of  power  spectrum  analysis  [2].  It  is  well  known  that  for  sinewaves  of  the 
same  strength  the  peak  depends  on  the  frequency  f:  highest  peak  for  f  =  mAf,  lowest 
for  f  =  (m+0.5)M,  where  m  is  an  integer.  This  means  that  a  scatterer  will  appear  larger 
if  it  is  in  the  center  of  the  FFT  cell  than  if  it  is  on  the  border  between  two  FFT  cells. 
Figure  5.3  shows  that  in  the  case  of  a  Kaiser  window,  the  peak  variation  is  about  25% 
which  is  an  unacceptable  ratio.  The  variation  in  peak  amplitude  due  to  the  FFT 
sampling  error  is  reduced  from  25%  to  5%  through  frequency  interpolation.  This  is 
accomplished  by  adding  N  zeros  to  the  N  return  samples  prior  to  FFT  transformation. 
The  effect  is  to  halve  Af.  The  number  of  significant  samples  without  interpolation  is  3 
and  with  interpolation  it  becomes  6.  We  have  shown  that  the  sum  of  the  significant 
samples  remains  practically  invariant,  with  or  without  interpolation. 

As  a  result  of  the  combined  proposed  corrections,  the  impulse  response  within 
one  footprint  becomes  invariant  within  ±4%.  This  is  observed  in  Fig.  12.1  which 
displays  the  reflectivity  of  a  footprint  versus  range  and  azimuth  for  eleven  scatterers 
distributed  over  the  footprint.  Therefore,  the  detection  of  scatterers  within  one  single 
footprint  can  now  be  considered  as  impulse  invariant. 

Invariant  mapping  is  used  to  transfer  the  footprint  information  onto  an  absolute 
x— y  ground  reflectivity  map.  In  this  manner  one  can  accumulate  the  information  for  a 
sequence  of  footprints  which  should  significantly  improve  the  detection  performance. 
The  mapping  technique  is  invariant  in  terms  of  volume.  In  other  words,  an  impulse  is 
mapped  onto  an  impulse  of  the  same  volume.  However,  the  impulse  obtained  after 
invariant  mapping  is  slightly  wider  and  shorter,  which  cannot  be  avoided.  The  peaks 
after  invariant  mapping  vary  within  ±5.5%  depending  on  the  relative  position  of  the 
FFT  cell  representing  the  scatterer  with  respect  to  the  rectangular  grid.  Figure  13.1 


20 


shows  that  the  x— y  reflectivity  map  obtained  using  invariant  mapping  is  nearly  impulse 
invariant 

7.  Other  Useful  Changes 

Several  major  improvements  to  the  RRDP  simulation  were  discussed  in  the 
previous  sections.  They  include:  (1)  extension  to  3— D  trajectories,  (2)  an  accurate 
circular  dish  simulation,  (3)  an  expanded  Taylor  expansion  which  includes  the  effect  of 
scatterer  height,  (4)  reduction  of  impulse  invariance  through  EFT  interpolation  and  a 
better  compensation  for  range  attenuation  and  antenna  gain.  This  section  discusses 
other  useful  changes,  except  for  those  related  to  graphics  displays  and  invariant 
mapping  which  are  documented  in  separate  sections. 

Two  main  programs  were  developed:  (1)  FRRDP,  where  the  output  is  a 
reflectivity  map  versus  range  and  azimuth  and  (2)  CRRDP,  where  the  output  is  a 
reflectivity  map  versus  absolute  x— y  coordinates.  Both  programs  have  the  same  input 
files:  CRRDPI.INP  which  defines  the  system  parameters,  and  CLUTINFI.INP  which 
defines  the  scatterers.  Also  subroutine  RECPOL  has  been  deleted.  The  overviews  of 
FRRDP  and  CRRDP  are  in  sections  9  and  11,  respectively. 

The  subroutine  CLUTINFO  which  is  used  to  enter  the  scatterer  data  was  also 
rewritten.  Major  changes  include:  (1)  the  position  of  C  can  be  entered  as  {x-.y-,^}  or 
{Hi^gi’V’  where  h,  is  the  scatterer  height,  (2)  the  output  array  SMOMAP  has  been 
modified  to  account  for  h.,  (3)  the  subroutines  XYTOMU  and  MUTOXY  were 
eliminated.  As  the  ith  scatterer  is  read  in  the  file,  the  information  {x-.y.,^^,^}  is 
displayed  on  the  screen.  CLUTINFO  provides  information  to  the  main  program,  on  the 
scatterers  that  are  both  within  the  selected  window  display  and  the  footprint,  through 
four  arrays:  SMOMAP,  WADRLO,  WADRUP,  and  WADSMO.  Range  smoothing 
performed  by  RANGSMOO,  transforms  one  real  scatterer  into  three  smoothed 
scatterers,  and  this  information  is  stored  in  SMOMAP  as  15  components: 
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{amp(l),  sfas(l),  mu(l),  jranon-1,  zsca(l),  amp(2),  sfas(2), 
mu(2),  jranon,  zsca(2),  amp(3),  sfas(3),  mu(3),  jranon+1,  zsca(3)}  (7.1) 

where  the  indexes  {1,2,3}  refers  to  {early,  on,  late},  and  jranon  is  the  range  index  for 
the  on-component.  The  elements  of  the  arrays  WADSMO  and  WJRAN  for  the  kth 
smoothed  scatterer  are:  {WADSMO(k)  =  address  of  the  scatterer  range  index  in 
SMOMAP}  and  { WJRAN (k)  =  SMOMAP(k)  =  range  index  for  scatterer}.  Once 
completed,  the  two  arrays  {WADSMO,  WJRAN}  are  sorted  with  respect  to  increasing 
range  index.  The  arrays  WADRLO  (JRAN)  and  WADRUP(JRAN),  where  JRAN 
varies  over  all  range  cells,  contain  the  lower  and  upper  addresses  of  the  sorted 
scatterers  which  are  within  the  JRAN  cell.  If  WADRLO  (JRAN)  equals  zero  this 
means  that  no  useful  scatterer  falls  within  this  ceil.  On  the  other  hand  {WADRLO 
(JRAN)  =  kl,  WADRUP  (JRAN)  =  k2}  means  that  the  sorted  scatterers  whose  index 
satisfies  kj  <  k  <  fall  within  the  JRAN  cell.  The  information  for  the  kth  sorted 
scatterer  can  be  found  in  SMOMAP  which  is  addressed  using  WADSMO(k). 

In  both  main  programs  (FRRDP  and  CRRDP)  and  in  CLUTINFO  the  inputs  are 
entered  interactively.  More  precisely,  the  steps  are  as  follows:  (1)  the  current  input 
data  file  is  read,  (2)  current  values  are  printed  on  the  screen,  (3)  a  prompt  K  =  ?  allows 
the  user  to  specify  the  index  of  the  data  to  be  changed,  (4)  if  the  3rd  data  input  need  to 
be  changed  to  200.0  we  enter  K  =  3,  then  200.0  at  the  X  =  ?  prompt,  and  (5)  when  all 
the  necessary  changes  have  been  made,  enter  K  =  0  and  the  updated  data  is  rewritten 
on  the  input  file.  This  approach  is  very  practical  because  usually  no  more  than  two 
parameters  are  changed  in  a  new  run.  Subroutines  WINPUT  and  WCHANGE  have 
been  developed  for  this  purpose,  they  are  listed  in  Appendix  E. 

Both  signal  synthesis  and  deramping  involve  very  large  phases.  These  phases 
have  to  be  evaluated  modulo  2tc.  We  found  that  it  is  critical  to  use  double  precision  for 
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phases  and  phase  residues.  This  affects  subroutines  TSIGNAL  and  DERAMP  which 
are  listed  in  Appendix  E. 

8.  Display  of  FRRDPG.GRD  and  CRRDPG.GRD  Using  SURFER 

The  main  programs  FRRDP  and  CRRDP  generate,  respectively,  the  files 
FRRDPG.GRD  and  CRRDPG.GRD.  The  data  in  FRRDPG.GRD  describes  the 
reflectivity  versus  azimuth  and  range.  The  data  in  CRRDPG.GRD  describes  the 
reflectivity  versus  absolute  x— y  coordinates.  The  purpose  of  this  section  is  to  explain 
how  the  format  of  FRRDPG.GRD  and  CRRDPG.GRD  was  selected  to  allow  direct 
display  by  SURFER  [7]. 

The  information  for  a  3— D  display  or  a  contour  map  could  be  written  as  a  three 

column  file,  where  the  number  of  rows  is  equal  to  the  number  of  data  points  and  one 

row  contains  the  x— y— z  coordinates  for  an  output  point.  The  disadvantage  of  such  a 

format  is  that  SURFER  will  require  gridding  prior  to  display.  Therefore,  it  is  much 

better  to  write  the  files  FRRDPG.GRD  and  CRRDPG.GRD  in  the  format  specified  by 

SURFER.  This  output  format  is  as  follows:  (1)  five  rows  which  specify  the  file  as 

indicated  by  the  graphics  manual,  and  (2)  the  data  in  blocks  of  length  N  ,  one  block  for 

each  constant  y.  There  are  a  total  of  blocks,  and  each  block  is  written  in  rows  of  10 

real  numbers.  Note  that  (N  ,N  )  are  the  dimension  of  the  x  and  y  coordinates. 

x  y 

SURFER  from  Golden  Software,  Inc.,  [7]  is  used  to  transform  CRRDPG.GRD 
into  a  three  dimensional  display.  In  a  standard  application  the  first  step  is  to  compute  a 
grid  starting  from  the  x— y— z  data.  However,  gridding  is  a  very  time  consuming 
process  and  it  can  cause  errors  if  the  smoothing  parameters  are  not  selected  properly. 
Fortunately,  when  the  increments  for  x  and  y  are  constant  with  no  data  missing,  the 
gridding  process  can  be  entirely  avoided  through  proper  formatting.  This  is  the  case  in 
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our  application.  It  follows  that  CRRDPG.GRD  is  written  in  the  block  format  specified 
in  the  previous  paragraph. 

The  file  CRRDPG.GRD  is  in  the  correct  format  for  a  3— D  display  using 
subroutine  SURF,  or  for  a  contour  map  using  subroutine  TOPO  [7].  Both  SURF  and 
TOPO  allow:  (1)  selection  of  color,  (2)  discrete  or  continuous  lines,  (3)  scaling,  (4) 
rotation,  (5)  labeling,  (6)  entering  a  title,  and  (7)  screen  gridding.  TOPO  also  provides 
control  of  the  minimum  and  maximum  contour  levels,  and  of  the  contour  increments. 
In  both  cases  we  can  obtain  a  hard  copy  by  using  one  of  four  devices:  (1)  an  IBM 
graphics  Proprinter,  (2)  an  HP  7475A  pen  plotter,  (3)  an  HP  laserjet  series  II  printer, 
and  (4)  an  HP  paint  jet  printer  which  allows  the  use  of  colors.  In  this  report  all  the 
hard  copy  displays  were  obtained  using  the  HP  laserjet  printer.  Examples  of  footprint 
reflectivity  maps  are  infection  13. 

9.  Overview  of  Computer  Simulation  for  Detection  Within  One  Footprint:  FRRDP 

The  overview  of  FRRDP  is  best  explained  using  the  functional  block  diagram  of 
Fig.  9.1.  The  goal  of  this  computer  simulation  is  to  write  the  file  FRRDPG.GRD  which 
SURFER  uses  to  display  the  reflectivity  map  versus  range  and  azimuth.  A  second 
output  file,  FRRDPP.PRT,  is  only  used  for  debugging.  Note  that  SURFER  includes 
both  SURF  and  TOPO.  The  block  diagram  is  an  implementation  of  the  techniques  and 
formulas  defined  in  the  previous  sections.  The  programs  are  written  in  Fortran  77,  in  a 
modular  structure  that  allows  easy  modifications.  They  are  listed  in  Appendix  E.  A 
functional  description  follows. 

The  initialization  block  includes  two  steps.  In  the  first  step  the  input  parameters 
are  entered  interactively  by  reading,  displaying  on  the  screen,  and  updating  the  input 
parameter  file  CRRDPI.INP.  A  printout  of  CRRDPI.INP  is  in  section  13.  Note  that 
the  input  data  includes  the  definition  of  the  display  window  through  the  array  WPAR. 
Only  the  scatterers  within  the  specified  window  are  involved  in  the  simulation  and 
mapping.  The  subroutine  WINPUT  is  used  to  enter  the  parameters  interactively  as 
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Figure  9.1  Block  diagram  of  the  FRRDP  simulation: 

SYNMAP  contains  reflectivity  versus  azimuth  and  range 


explained  in  section  7,  and  PRTXY  displays  the  current  parameters  on  the  screen. 
Subroutine  WINPUT  calls  WCHANGE  which  uses  KINPUT  and  XINPUT  to  enter 
integer  or  real  numbers  on  line.  The  second  step  is  the  precomputation  of  the  trajectory 
and  footprint  information  at  central  dwell  time.  Also  the  output  reflectivity  array 
SYNMAP  is  initialized.  The  array  WMU3DB  stores  the  maximum  azimuth  increments 
pmx,  one  for  each  range  cell. 

The  ground  map  is  computed  using  subroutine  CLUTINFO.  The  scatterer  data 
is  in  the  file  CLUTINFI.INP.  The  format  is  as  follows:  (1)  the  first  line  is  the  signed 
number  of  scatterers,  where  a  positive  sign  indicates  x— y  coordinates  for  the  scatterer, 
while  a  negative  sign  means  relative  spherical  coordinates,  and  (2)  each  scatterer  is 
entered  as  three  lines  of  data,  either  (x.,y  h,)  or  (u.,r-,h-).  The  scatterer  data  is  read, 
displayed  on  the  screen,  modified  on  line,  and  rewritten  on  the  file  CLUTINFI.INP, 
using  subroutine  WINPUT  which  was  discussed  earlier.  Only  the  scatterers  that  fall 
within  the  selected  footprint  and  display  map  (including  range  smoothing  effect)  are 
processed  further.  As  well  known,  a  point  scatterer  affects  three  consecutive  range 
cells  because  of  the  finite  bandwidth  of  the  system.  This  effect  which  is  denoted  as 
range  smoothing  is  simulated  using  subroutine  RANGSMOO,  which  splits  each 
scatterer  into  three  equivalent  scatterers.  All  the  necessary  scatterer  information  is 
passed  to  the  main  program  through  the  arrays  SMOMAP,  WADRLO,  WADRUP,  and 
WADSMO.  More  details  on  CLUTINFO  are  in  section  7.  It  is  convenient  to  use  the 
index  KNPPP  which  varies  from  1  to  NFFT  from  minimum  to  maximum  doppler, 
rather  than  the  FFT  index  KFF  which  does  not  correspond  to  a  continuous  doppler 
variation.  The  correspondence  is:  {KFF  =  <KNPPP  +  NFFT/2  -  1  >NFFT  +  1,  KNPPP 
=  1,  NFFT},  where  <  >NFFT  means  modulo  NFFT.  Furthermore,  the  useful  range  of 
KNPPP  is  limited  to  within  the  3dB  footprint,  that  is  (KPPPMN  <  KNPPP  < 
NPPPMX}.  The  change  from  KNPPP  to  KFF  is  effected  by  the  array  LNPPP. 
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The  heart  of  FRRDP  is  loop  1  which  corresponds  to  "DO  35"  in  the  program. 
In  this  loop  the  range  cells  are  processed  sequentially:  JRAN  =  1,  JRANMX. 

In  Block  3  a  check  for  the  presence  of  smoothed  scatterers  within  the  current 
JRANth  cell  is  performed.  If  WADRLO(JRAN)  =  0  there  are  no  scatterers  and  loopl 
has  been  completed. 

Block  4  represents  loop  2  which  corresponds  to  "DO  30"  in  the  program.  Its 
purpose  is  to  synthesize  the  signal  (XD(KDQ),  XQ(KDQ),  KDQ  =  1,  NFFT}  in  the 
JRANth  range  cell.  The  contributing  scatterers  are  defined  by  { WADRLO(JRAN)  < 
KTHSCA  <  WADRUP(JRAN) } .  The  information  for  the  KTHSCA  scatterer  is 
transferred  from  SMOMAP  into  the  array  WTSMO.  The  antenna  gain  and  the 
contribution  to  the  signal  due  to  the  KTHSCA  scatterer  are  computed  using  subroutines 
SCAGAN  and  TSIGNAL  which  are  listed  in  Appendix  E. 

The  electronic  processing  of  the  JRANth  cell  is  in  block  5.  It  consists  of  the 
following  steps:  (1)  deramping,  (2)  windowing,  and  (3)  fast  Fourier  transformation 
(FFT).  The  subroutines  DERAMP,  WINDOWS  and  FFT  which  are  listed  in 
Appendix  E  are  used  for  this  purpose.  Note  that  the  FFT  is  performed  on  the  extended 
arrays  (XD(k),  XQ(k),  k  =  1,  2*NFFT)  which  are  padded  with  NFFT  zeros.  This 
corresponds  to  the  implementation  of  frequency  interpolation  which  was  discussed  in 
section  6. 

Block  6  is  concerned  with  scaling  of  the  FFT  cells  and  the  synthesis  of  the 
reflectivity  array  SYNMAP.  First,  the  array  WAGNCL  which  contains  the  antenna 
gains  at  the  center  of  the  FFT  cells  is  computed,  ('T)0  39"  in  the  program).  Then,  both 
the  scaling  and  the  computation  of  SYNMAP  are  performed  in  "DO  41".  The  scaling 
which  compensates  for  antenna  gain  and  range  variation  is  the  implementation  of 
section  6. 
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Block  7  coiresponds  to  the  two  "DO  60"  in  the  program.  The  purpose  of  this 
block  is  to  write  the  file  FRRDPG.GRD  which  contains  the  reflectivity  map  in  the 
format  suggested  by  [7]  (see  section  8). 

10.  New  Invariant  Mapping  Technique 

Application  of  SURF  on  the  file  FRRDPG.GRD  written  by  FRRDP  produces  a 
3— D  reflectivity  map,  where  {x  =  azimuth,  y  =  range,  z  =  reflectivity}.  In  many 
applications  it  is  desirable  to  represent  the  information  as  an  absolute  reflectivity  map 
[x=  absolute— x,  y  =  absolute— y,  z  =  reflectivity}.  This  approach  allows  combining  the 
information  of  a  sequence  of  footprints  into  an  overall  absolute  map.  An  invariant 
mapping  technique  to  transfer  the  information  from  the  footprint  to  an  absolute  map 
was  developed  in  [1,2].  The  reflectivity  of  one  trapezoidal  FFT  cell  is  shared  between 
the  overlapping  rectangular  grid  cells  on  the  basis  of  overlap  area.  It  follows  that  the 
efficiency  of  the  invariant  mapping  algorithm  depends  mostly  on  the  computation  of 
the  area  of  overlap.  In  [1,2]  the  trapezoid  was  split  into  two  triangles  and  various 
overlap  configurations  were  considered.  In  this  report  we  are  presenting  a  new 
algorithm  which  computes  directly  the  overlap  of  the  trapezoid  cells  with  the  absolute 
grid  rectangles.  The  new  algorithm  is  much  more  general,  compact,  and  efficient  than 
that  of  [1,2].  It  is  discussed  below. 

Prior  to  writing  the  algorithm  we  need  to  define:  (1)  the  trapezoidal  FFT  cells 
being  mapped  and  their  reflectivity,  (2)  the  absolute  rectangular  cells  and  their 
estimated  reflectivity,  and  (3)  the  area  of  overlap  between  the  mth  trapezoidal  cell  and 
the  kth  rectangle; 

TRP(m)  <=>  mth  trapezoidal  cell 

RFT(m)  <=>  reflectivity  of  the  mth  trapezoidal  cell 

REC(k)  <=>  kth  absolute  rectangular  cell 

RFR(k)  estimated  reflectivity  of  the  kth  rectangular  cell 

OVL(m,k)  <=>  overlap  area  between  mth  trapezoid  and  kth  rectangle. 
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Using  the  above  definitions,  the  invariant  mapping  algorithm  implements  the 
following  formula 


IOVL(ra,k)  RFT(m) 

RFR(k)  =  — m -  (10.1) 

lOVL(m,k) 


where  m  defines  all  the  trapezoids  which  overlap  the  kth  rectangle. 

The  invariant  mapping  algorithm  is  implemented  in  CRRDP  which  is  discussed 
in  section  11.  It  involves  the  following  subroutines:  APEXTRAP  and  MAPRGM.  The 
subroutine  MAPRGM  calls  MAPAZM  which  calls  POLONFRA.  The  only  complicated 
subroutine  is  POLONFRA  which  is  discussed  next. 

Figure  10.1  shows  the  mth  trapezoid  being  mapped  on  its  frame  prior  to  transfer 
to  the  absolute  window  defined  by  WPAR  =  { DELMAP.THRESHXMN,  YMN,XMX, 
YMX}.  Only  the  frame  that  contains  the  trapezoid  is  involved  in  the  overlap.  The 
purpose  of  POLONFRA  is  to  compute  all  the  nonzero  OVL(m,k).  The  inputs  of 
POLONFRA  are  as  follows: 

NBA  *=>  #  trapezoid  apexes  =  4 

WAX  <=*  x— coordinates  of  trapezoid  apexes 

WAY  £3  y— coordinates  of  trapezoid  apexes 

WPAR  «  array  defining  the  absolute  window  being  mapped. 

The  outputs  of  POLONFRA  are: 

(HOJJO)  <=*  absolute  coordinates  of  frame  upper  left  comer 
(IMXJMX)  «=>  frame  dimensions 
TAREA  area  of  trapezoid 

WWAREA(I,J)  «  overlapping  area  for  the  (I,J)th  rectangle  within  the  frame. 
IOKRO  =  0  if  apexes  are  entered  in  counterclockwise  order,  and  1  otherwise. 
EFAEL  $=*  test  for  failure,  it  should  be  equal  to  zero. 
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(XMN.YMN)-upper  left  corner  of  window 
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Figure  10.1  Mapping  the  mth  trapezoidal  cell  on  its  frame 
prior  to  transfer  to  absolute  grid  map. 
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Conceptually  POLONFRA  operates  as  follows:  (1)  following  the  trapezoid 
counterclockwise  contour  (CCWC),  one  computes  all  the  intersects  with  the  rectangular 
grid,  (2)  the  intersects  are  ordered  in  order  of  CCWC,  (3)  duplicate  intersects  are 
eliminated,  (4)  the  comers  of  the  rectangles  within  the  frame  are  defined  as  exterior  or 
interior  to  the  trapezoid,  and  (5)  the  overlap  areas  of  the  rectangles  which  have  interior 
and/or  intersects  with  the  trapezoid  are  computed.  The  overlap  contour  for  the 
rectangle  is  defined  by  ordering  the  interior  comers  and/or  intersect  points  in  CCWC. 
In  this  manner  we  define  a  convex  contour  whose  area  is  easy  to  compute  using  the 
function  AREPOL.  Listings  for  POLONFRA  and  associated  subprograms  are  in 
Appendix  F. 

11.  Overview  of  Computer  Simulation  with  Invariant  Mapping:  CRRDP 

The  overview  of  CRRDP  is  best  explained  using  the  functional  block  diagram 
of  Fig.  11.1.  The  goal  of  this  computer  simulation  is  to  write  the  file  CRRDPG.GRD 
which  SURFER  uses  to  display  the  reflectivity  map  versus  the  absolute  x— y 
coordinates.  A  second  output  file,  CRRDPP.PRT,  is  only  used  for  debugging.  The 
main  difference  between  FRRDP  and  CRRDP  is  that  FRRDP  generates  the  footprint 
information  directly,  while  CRRDP  transfer  the  information  onto  a  absolute  x— y  grid 
using  the  invariant  mapping  technique  discussed  in  section  10. 

It  follows  that  the  first  block  of  CRRDP  is  the  same  as  the  first  block  of  FRRDP 
except  for  the  initialization  of  SYNMAP  and  SUMARA  ("DO  13”  is  the  program). 
These  two  arrays  are  used  in  subroutine  MAPRGM  to  accumulate  both  the  numerator 
and  denominator  of  (10.1).  Blocks  {2,3,4, 5}  in  CRRDP  are  the  same  as  in  FRRDP  and 
were  discussed  in  section  9. 

Block  6  performs  antenna  and  range  attenuation  compensation  as  explained  in 
section  9.  Block  7  implements  Invariant  Mapping.  It  performs  two  tasks:  (1) 
computation  of  the  trapezoid  FFT  cells  apexes  coordinates  in  the  current  range  cell 
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Figure  11.1  Block  Diagram  for  the  CRRDP  Simulation: 

SYNMAF  contains  the  absolute  x-y  reflectivity  map. 
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(subroutine  APEXTRAP),  and  (2)  computation  of  the  arrays  SYNMAP  and  SUMARA 
which  represent  respectively,  the  numerator  and  denominator  of  (10.1).  The 
accumulation  in  SYNMAP  and  SUMARA  is  performed  by  subroutine  MAPAZM  which 
is  called  by  MAPRGM. 

Block  8  corresponds  to  the  two  "DO  60"  in  the  program.  First,  the  reflectivity 
map  SYNMAP  is  computed  according  to  (10.1).  Then,  the  reflectivity  map  versus 
absolute  x— y  coordinates  is  written  on  the  file  CRRDPG.GRD  using  the  format 
recommended  by  [7],  Appendix  G  contains  the  listings  for  CRRDP  and  its  subroutines, 
except  for  those  listed  earlier  in  Appendixes  C,  D,  E,  and  F. 

12.  Demonstration  of  the  FRRDP  Simulation 

In  order  to  run  FRRDP  or  CRRDP  one  must  first  define  the  following:  (1)  the 
seeker  antenna  information,  (2)  the  footprint  center,  (3)  the  selected  frame  window,  and 
(4)  the  electronic  processing  parameters.  This  is  accomplished  interactively  using  the 
subroutine  WINPUT  on  the  data  file  CRRDPI.INP.  Table  12.1  shows  the  file 
CRRDPI.INP  with  data  set  #1,  which  is  used  for  most  of  our  examples.  This  file  is  self 
explanatory  since  each  parameter  is  followed  by  a  comment.  A  demonstration  diskette 
(DEMO  #1)  is  included  in  [8]. 

The  clutter  configuration  is  also  entered  interactively  through  subroutine 
CLUTINFO  which  uses  WINPUT  to  read,  and  update  as  needed,  the  clutter  file 
CLUTINR.INP.  This  file  contains  the  signed  number  of  scatterers  followed  by  (p-Pg., 
hj)  for  negative  or  {x-,  y-,  h- }  for  positive  sign  (see  section  7).  Table  12.2a  contains 
the  data  for  eleven  scatterers  with  zero  height  which  are  distributed  over  the  footprint 
specified  by  Table  12.1.  When  FRRDP  is  executed  using  Table  12.2a  the  scatterer 
information  appears  on  the  screen  as  shown  in  Table  12.2b. 

All  the  3— D  displays  in  this  section  are  for  a  312°  rotation  about  the  z— axis 
(entered  in  SURF).  Figure  12.1  shows  the  3— D  reflectivity  map  versus  azimuth  and 
range,  obtained  by  using  SURF  [7]  on  the  file  FRRDPG.GRD  created  by  FRRDP  for 
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Table  12.1.  Data  set  #1  for  file  CRRDPI.INP 


.0100000 
.0200000 
150.0000000 
.0000000 
1.2217300 
.0000000 
7.0000000 
900.0000000 
516.2463000 
744.2216000 
.0000000 
.0000000 
.0000000 
. 0000000 
1.0000000 
6.0000000 
.7500000 
. 1000000 
480.0000000 
710.0000000 
550.0000000 
780.0000000 
.0000000 
4.0000000 
3.1415900 
2.0000000 
1.0000000 


TIMCNT=central  time 

TIMDWL=observation  interval 

VELOM=seeker  velocity 

ALFM=seeker  azimuth  angle 

BETM=seeker  elevation  angle 

WMPOSC(l)=seeker  x-coordinate 

WMPOSC (2) =seeker  y-coordinate 

WMPOSC (3) “seeker  z-coordinate 

WQ9 (1) “footprint  center  x-coordinate 

WQ9 (2) “footprint  center  y-coordinate 

DUMMY=  future  use 

DUMMY=future  use 

DUMMY= future  use 

DUMMY= future  use 

DELRAN= range  bin 

FFNBIT=#  FFT  bits 

WPAR(l) “rectangular  grid  increment 
WPAP (2) “mapping  threshold 
WPAR ( 3 ) =x-min  display  window 
WPAR ( 4 ) =y-min  display  window 
WPAR ( 5 ) =x-max  display  window 
WPAR (6) =y-max  display  window 
XIFXMT“transmitter  phase 
SELWIN“define  FFT  window 
BETKAI=Kaiser  window  parameter 
SELDIS=*define  dish  simulation 
BITINT“define  FFT  interpolation 
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Table  12.2a.  Data  set  #1  for  input  file  CLUTINFI . INP . 

-34 . 00000000  Wl=3*  #scatt . +1 , pos .  If  (xi,yi) , neg.  If  (mui , rgi) 

.0000000  Xl=x(lst  scat),  or  mul 

1272.7910000  Yl=y(lst  scat),  or  rgl 

0.0000000  Zl3z ( 1st  scat)  *  hbarl 
-.0304062  Xl*x(2nd  scat),  or  mu2 - 

1272.7910000  Yl=y(2nd  scat),  or  rg2 

.0000000  Zl=«z(2nd  scat)  3  hbar2  • 

.0304062  Xl*»x(3rd  scat),  or  mu3  • 

1272.7910000  Yl*y(3rd  scat),  or  rg3 

.0000000  Zl»z(3rd  scat)  *  hbar3 
.0000000  Xl3x(4th  scat),  or  mu4 

1285.7910000  Yl*y(4th  scat) ,  or  rg4 

.0000000  Zl=z(4th  scat)  ®  hbar4 

-.0267718  Xl=x(5th  scat),  or  mu5 

1285.7910000  Yl*y(5th  scat),  or  rg5 

.0000000  Zl=z(5th  scat)  =®hbar5 
.0267718  Xl=x(6th  scat),  or  mu6 

1285.7910000  Yl=y(6th  scat) ,  or  rg6 

.0000000  Zl=z(6th  scat)  ®  hbar6 
.ooooooo  Xl=x(7th  scat),  or  mu7 

1260.7910000  Yl=y (7th  scat),  or  rg7 

.0000000  Zl3z(7th  scat)  3  hbar7 

-.0273642  Xl3x(8th  scat),  or  mu8 

1260.7910000  Yl3y(8th  scat),  or  rg8 

.0000000  Zl3z(8th  scat)  3  hbar8 

.0273642  Xl=x ( 9th  scat),  or  mu9 

1260.7910000  Yl=y(9th  scat) ,  or  rg9 

.ooooooo  Zl3z(9th  scat)  3  hbar9 

.0000000  Xl3x(10th  scat),  or  muio 

1295.7910000  Yl=y(10th  scat) ,  or  rglO 

.ooooooo  Zl3z(l0th  scat)  3  hbarlO 
.0000000  Xl3x(llth  scat),  or  mull 

1250.7910000  Yl=y(llth  scat) ,  or  rgll 

.0000000  Zl3z(llth  scat)  3  hbarll 
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Table  12.2b.  Screen  printout  for  scatterer  input  in  Table  12.2a 


XI* . 51624E+03 
XI* . 49359E+03 
XIs®.  53842E+03 
XI*. 52674E+03 
XI* . 50641E+03 
XI3®.  54668E+03 
XI*.50646E+03 
XI* . 48648E+03 
XIs®.  52606E+03 
XIs®.  53474E+03 
XI*. 49824E+03 


YI* . 74422E+03 
YI*. 75957E+03 
YI*. 72818E+03 
YI*. 75920E+03 
YI*. 77303E+03 
YI*. 74483E+03 
YI*. 73025E+03 
YI3® .  74384E+03 
YI*.  71612E+03 
YI3®.  77063E+03 
YI®* .  71850E+03 


MUI= 
MUI= 
MU  Is* 
MUI= 
MUI* 
MUI®® 
MUI* 
MUI* 
MUI* 
MUI* 
MUI* 


. 00000E+00  RGI*. 12728E+04  HI* 
- . 30406E-01  RGI*. 12728E+04  HI* 
. 30406E-01  RGI*. 12728E+04  HI* 
. 00000E+00  RGI*. 12858E+04  HI* 
-. 26772E-Q1  RGI*. 12858E+04  HI* 
. 26772E-01  RGI*. 12858E+04  HI* 
. OOOOOE+OO  RGI*. 12608E+04  HI* 
-. 27364E-01  RGI*. 12608E+04  HI* 
.27364E-01  RGI*. 12608E+04  HI* 
.00000E+00  RGI*. 12958E+04  HI* 
.00000E+00  RGI*. 12508E+04  HI* 


. 00E+00 
.  00E+00 
.  00E+00 
.  00E+00 
.  00E+00 
.  00E+00 
.  00E+00 
.  00E+00 
.  00E+00 
.  00E+00 
.  OOE+OO 
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the  inputs  defined  in  Tables  12.1  and  12.2a.  Figure  12.2  shows  a  side  view  of  the  3— D 
display  in  Fig.  12.1.  This  demonstrates  that  impulse  invariance  has  been  achieved 
within  ±4%  over  the  entire  footprint,  and  it  indicates  that  antenna  and  range 
compensation  work  very  well.  Figure  12.3  is  a  reflectivity  contour  plot  obtained  by 
using  TOPO  [7]  on  the  same  file  FRRDPG.GRD.  This  figure  shows  that:  (1)  the 
scatterers  are  detected  at  the  correct  locations,  and  (2)  the  azimuth  resolution  is  now 
comparable  to  range  resolution. 

Even  though  the  footprint  defined  by  Table  12.1  contains  eleven  scatterers 
defined  by  Tables  12.2a  and  12.2b,  the  reflectivity  display  can  be  restricted  to  a  single 
scatterer  by  shrinking  the  size  of  the  display  window  (WPAR)  in  CRRDPi.iN^.  Table 
12.3  is  identical  to  Table  12.1  except  for  a  smaller  window  about  the  center  of  the 
footprint.  The  output  reflectivity  for  the  run  that  uses  Tables  12.3  and  12.2a  is  shown 
in  Fig.  12.4.  To  demonstrate  the  effect  of  scatterer  height,  =  3  meters  is  entered  for 
the  central  scatterer  as  shown  in  Table  12.4.  The  reflectivity  displays  for  the  run  using 
Tables  12.3  and  12.4  are  in  Figs.  12.5a  and  b.  This  experiment  demonstrates:  (1)  the 
utilization  of  the  window  display,  and  (2)  the  effect  of  scatterer  height.  Indeed  a 
scatterer  height  of  3m  has  the  effect  of  reducing  range  and  increasing  doppler  by  about 
one  cell. 

13.  Demonstration  of  CRRDP  Simulation 

The  examples  in  this  section  use  the  same  inputs  as  those  defined  for  FRRDP  in 
the  previous  section.  As  before,  the  system  and  the  clutter  parameters  are  in  files 
CRRDPI.INP  and  CLLl  llNFI.INP,  respectively,  and  are  entered  interactively  using 
WTNPUT.  The  examples  in  this  section  are  similar  to  those  in  the  previous  section 
except  that  the  outputs  are  absolute  reflectivity  maps  versus  absolute  x-  y  coordinates. 
A  demonstration  diskette  (DEMO  #2)  is  included  in  [8] 

All  the  3— D  displays  in  this  section  are  for  a  312°  rotation  about  the  z— axis 
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1.00  12  00  23  00  34  00  45  00  50  00  67  00  70  00  69  00  100.00  11100  122  00 


Table  12.3.  Data  set  #2  for  file  CRRDPI.INP 


.0100000 
. 0200000 
150.0000000 
.0000000 
1.2217300 
.0000000 
7.0000000 
900.0000000 
516.2463000 
744.2216000 
.0000000 
.0000000 
.0000000 
.0000000 
1.0000000 
6.0000000 
.7500000 
. 1000000 
506.0000000 
734.0000000 
526.0000000 
754.0000000 
.0000000 
4.0000000 
3.1415900 
2.0000000 
1.0000000 


TIMCNT=central  time 

TIMDWL=observation  interval 

VELOM=seeker  velocity 

ALFM=seeker  azimuth  angle 

BETM=seeker  elevation  angle 

WMPOSC(l)=seeker  x-coordinate 

WMPOSC (2) =seeker  y-coordinate 

WMPOSC ( 3 ) =seeker  z-coordinate 

WQ9 (1) =footprint  center  x-coordinate 

WQ9 (2) “footprint  center  y-coordinate 

DUMMY® future  use 

DUMMY® future  use 

DUMMY® future  use 

DUMMY® future  use 

DELRAN=range  bin 

FFNBIT®#  FFT  bits 

WPAR(l) =rectangular  grid  increment 
WPAR ( 2 ) =mapping  threshold 
WPAR ( 3 ) =x-min  display  window 
WPAR ( 4 ) =y-min  display  window 
WPAR ( 5 ) =x-max  display  window 
WPAR ( 6 ) =y-max  display  window 
XIFXMT=transmitter  phase 
SELWIN=def ine  FFT  window 
BETKAI=Kaiser  window  parameter 
SELDIS=define  dish  simulation 
BITINT=define  FFT  interpolation 
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Table  12.4.  Data  set  #2  for  input  file  CLUTINFI . INP . 

-34 . 00000000  Wl=3*  #scatt.+l,pos.  If  (xi,yi) ,neg.  if  (mui ~ rgi) 

.0000000  Xl=x(lst  scat),  or  mul 

1272.7910000  Yl=y(lst  scat),  or  rgl 

3.0000000  Zl=z(lst  scat)  =  hbarl 
-.0304062  Xl=x(2nd  scat),  or  mu2 

1272.7910000  Yl=y(2nd  scat),  or  rg2 

.0000000  Zl=z(2nd  scat)  =  hbar2 
.0304062  Xl=x(3rd  scat),  or  mu3 

1272.7910000  Yl=y(3rd  scat),  or  rg3 

.0000000  Zl=z(3rd  scat)  =  hbar3 
.0000000  Xl=x(4th  scat),  or  mu4 

1285.7910000  Yl=y(4th  scat),  or  rg4 

.0000000  Zl=z(4th  scat)  =  hbar4 

-.0267718  Xl=x(5th  scat),  or  mu 5 

1285.7910000  Yl=y(5th  scat),  or  rg5 

.0000000  Zl=z(5th  scat)  =hbar5 
.0267718  Xl=x(6th  scat),  or  mu 6 

1285.7910000  Yl=y(6th  scat),  or  rg6 

.ooooooo  Zl=z (6th  scat)  =  hbar6 
.ooooooo  Xl=x(7th  scat),  or  mu7 

1260.7910000  Yl=y (7th  scat),  or  rg7 

.ooooooo  Zl=z (7th  scat)  =  hbar7 

-.0273642  Xl=x(8th  scat),  or  mu8 

1260.7910000  Yl=y(8th  scat),  or  rg8 

.0000000  Zl=z(8th  scat)  =  hbar8 
.0273642  Xl=x(9th  scat),  or  mu9 

1260.7910000  Yl=y(9th  scat),  or  rg9 

.ooooooo  zi=z(9th  scat)  =  hbar9 

.ooooooo  xi=x(loth  scat),  or  muio 

1295.7910000  Yl=y ( loth  scat) ,  orrglO 

.ooooooo  Zl=z(l0th  scat)  =  hbario 
.0000000  Xl=x(llth  scat),  or  mull 

1250.7910000  Yl=y(llth  scat),  or  rgll 

.ooooooo  Zl=z(llth  scat)  =  hbarll 
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(entered  in  SURF).  Figures  13.1a  and  b  shows  the  3-D  reflectivity  and  the  contour 
map  versus  absolute  x— y  coordinates  obtained  by  using  SURFER  [7]  on  the  file 
CRRDPG.GRD  created  by  CRRDP  for  the  inputs  defined  in  Tables  12.1  and  12.2a. 
Figure  13.2  shows  a  side  view  of  the  3— D  plot  in  Fig.  13.1a.  As  expected  the  invariant 
mapping  technique  produces:  (1)  a  small  widening  of  the  pulse,  and  (2)  a  slight 
increase  in  the  variation  over  the  footprint;  this  is  due  to  different  configurations  of  the 
FFT  trapezoidal  cells  relative  to  the  rectangular  grid.  However,  one  can  clearly 
observe  that:  (1)  the  peak  reflectivity  is  nearly  invariant,  (2)  the  maximum  variation 
over  the  x— y  map  is  within  ±5.5%,  which  should  be  acceptable  for  most  applications, 
and  (3)  the  volume  invariance  of  the  impulse  response  is  nearly  perfect;  this  provides 
some  compensation  during  detection  because  the  shorter  impulses  are  wider  than  the 
taller  impulses. 

As  earlier  one  can  reduce  the  selected  window  display  to  a  single  scatterer 
output  by  shrinking  the  display  window  (WPAR).  Figure  13.3  is  for  the  inputs 
specified  by  Tables  12.2a  and  12.3.  The  effect  of  scatterer  height  is  demonstrated  by 
running  CRRDP  with  the  inputs  defined  by  Tables  12.2a  and  12.4. 

All  the  results  so  far  were  for  rectangular  grids  with  grid  size 
(Ax  =  Ay  =  0.75m).  The  performance  of  the  invariant  mapping  technique  becomes 
even  better  when  the  grid  increment  is  reduced  to  0.6m.  Figure  13.4  shows  the 
reflectivity  for  the  input  files  given  in  Tables  13.1  and  12.2a.  Only  eight  scatterers  are 
showing  because  we  were  forced  to  decrease  the  size  of  the  of  the  display  window  due 
to  the  restriction  of  SURFER  to  a  maximum  of  10000  points.  The  impulses  are 
smoother  and  the  peak  variation  over  the  display  is  now  only  ±4%.  Conversely,  if  the 
size  of  the  grid  is  increased  to  1m,  the  impulse  variation  will  increase  to  about  ±10%. 
Note  that  in  all  cases  the  average  peak  value  remain  constant,  which  confirms  the 
invariance  of  the  mapping  technique. 
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from  CRRDP  (  Tab  I  as  12.1  &12.2). 


63  00  70  00  77  00  84 
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Table  13.1.  Data  set  #3  for  file  CRRDPI.INP 


.0100000 
.0200000 
150.0000000 
.0000000 
1.2217300 
.0000000 
7.0000000 
900.0000000 
516.2463000 
744.2216000 
.0000000 
. 0000000 
.0000000 
.0000000 
1.0000000 
6.0000000 
.7500000 
. 1000000 
491.0000000 
721.0000000 
550.0000000 
780.0000000 
.0000000 
4.0000000 
3.1415900 
2.0000000 
1.0000000 


TIMCNT=central  time 

TIMDWL=observation  interval 

VELOM= seeker  velocity 

ALFM=*seeker  azimuth  angle 

BETM=seeker  elevation  angle 

WMPOSC ( 1) =seeker  x-coordinate 

WMPOSC ( 2 ) =seeker  y-coordinate 

WMPOSC ( 3 ) =seeker  z-coordinate 

WQ9 (1) =footprint  center  x-coordinate 

WQ9 ( 2 ) =f ootprint  center  y-coordinate 

DUMMY=future  use 

DUMMY= future  use 

DUMMY= future  use 

DUMMY= future  use 

DELRAN=range  bin 


FFNBIT=#  FFT  bits 

WPAR (1) =rectangular  grid  increment 
WPAR ( 2 ) =mapping  threshold 
WPAR ( 3 ) =x-min  display  window 
WPAR (4) =y-min  display  window 
WPAR ( 5 ) =x-max  display  window 
WPAR(6)=y-max  display  window 
XIFXMT=transmitter  phase 
SELWlN=def ine  FFT  window 
BETKAI=Kaiser  window  parameter 
SELDIS=define  dish  simulation 
BITINT=define  FFT  interpolation 
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14.  Summary,  Conclusions  and  Recommendation  for  Future  Work 


This  report  extends  and  upgrades  the  computer  simulation  for  Range  Relative 
Doppler  Processing  (RRDP)  and  invariant  mapping  for  MM  Wave  Seekers.  The  main 
program  and  most  of  the  subroutines  are  significantly  different  from  those  in  [1,2]. 
Major  improvements  include:  (1)  a  more  general  geometry  where  the  trajectory  is  not 
restricted  to  a  vertical  plane,  (2)  an  exact  simulation  for  the  dish  antenna,  (3)  a  more 
exact  signal  synthesis  that  includes  the  effect  of  scatterer  height,  (4)  the  development  of 
an  interactive  input  subroutine  which  reads,  displays  on  screen,  prompts  for  online 
changes,  and  updates  the  data  file,  (5)  a  much  better  compensation  for  antenna  gain  and 
range  attenuation,  (5)  a  reduction  of  the  peak  reflectivity  variation  due  to  FFT 
quantization  through  frequency  interpolation,  (7)  a  new,  better,  and  more  efficient 
invariant  mapping  algorithm,  and  (8)  online  3— D  and  contour  graphic  displays  for  both 
types  of  reflectivity  maps. 

Two  main  programs  were  developed,  FRRDP  and  CRRDP.  The  first,  FRRDP, 
generates  the  file  FRRDPG.GRD  which  SURFER  uses  to  display  the  footprint 
reflectivity  map  versus  azimuth  and  range.  This  program  should  be  useful  when  one 
needs  to  perform  target  detection  within  one  footprint.  The  second  program,  CRRDP, 
generates  the  file  CRRDPG.GRD  which  SURFER  uses  to  display  the  reflectivity  versus 
the  absolute  x— y  grid.  Our  invariant  mapping  technique  is  used  to  transfer  the 
information  from  the  footprint  to  the  absolute  map.  This  program  should  be  useful  to 
accumulate  the  information  of  a  sequence  of  footprints  prior  to  final  detection. 

Examples  show  various  footprint  reflectivity  maps  obtained  by  using  FRRDP. 
These  examples  demonstrate  the  following:  (1)  the  scatterers  are  detected  at  the  proper 
locations,  (2)  impulse  invariance  has  been  achieved  over  the  entire  footprint,  and  (3) 
the  effect  of  scatterer  height  is  significant  and  has  been  properly  simulated. 
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Examples  also  show  various  absolute  x— y  reflectivity  maps  obtained  by  using 
CRRDP.  They  demonstrate  the  performance  of  Invariant  Mapping:  (1)  the  peaks 
remain  practically  the  same  in  value  and  location,  (2)  impulse  invariance  has  been 
achieved  over  the  absolute  map  within  a  ±5.5%  variation,  and  (3)  the  volume 
invariance  of  the  impulses  is  nearly  perfect.  Most  of  the  examples  were  run  with  a 
rectangular  grid  of  size  Ax  =  Ay  =  0.75  m.  Increasing  the  grid  size  to  lm  will  increase 
the  peak  variation  over  the  absolute  map  to  about  ±10%.  On  the  other  hand,  a 
reduction  of  the  grid  size  to  0.6m  reduces  the  variation  to  about  ±4%,  but  it  increases 
the  computation  requirements.  At  the  present  time  the  invariant  mapping  involves  only 
a  single  footprint  but  it  will  be  easy  to  extend  it  to  a  sequence  of  footprints.  One 
expects  that  the  impulse  invariance  over  the  absolute  map  would  become  even  better 
due  to  the  smoothing  effect  of  footprint  integration. 

Recommendations  for  future  work  include  the  following:  (1)  using  realistic 
target  models  as  input,  (2)  using  realistic  clutter  models  with  or  without  targets,  (3) 
incorporating  pulse  compression  in  the  simulation  for  a  significant  increase  in  the 
signal-to-noise  ratio,  (4)  adding  detection  and  pattern  recognition  to  the  programs,  and 
(5)  extending  the  simulation  to  multiple  polarization. 
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APPENDIX  A 

Relevant  Footprint  Information  For  m  =  {dx>  d^,  h) 

•»4  ■+ 

A.1  Footprint  center  q  =  q(t  ),  see  Fig.  A.l. 

V 

Rectangular  coordinates  of  footprint  center  =  {x^,  y^,  0} 
range:  r  *  (h2  +  (xq  -  dx)2  +  (yq  -  dy)2) 

*  i 

azimuth  angle:  a  =  tan  [(x  -d  )/(y  -d  )] 

4  a  4  / 

elevation  angle  =  (3  =  cos-1  — ^ 

*■  r  • 

A.2  Range  bounds  and  indexing,  see  Fig.  A.l. 

Minimum  range:  Ranmn  = - j - 

cos(p  -  8h) 

Maximum  range:  Ranmx  = - ^ - 

cos(J3  +  0H) 

0H  =  half  beam  width 

-♦*  *  * 

Range  index  for  q  :  j  =  INT[(r  —  Ranmn  +  Ar/2)/Ar] 

* 

Number  of  range  cells:  2j 

A.3  Information  for  the  jth  range  cell 

*  *  * 

median  range:  r  =  r  +  (j  -  j  )  Ar 
elevation  angle:  (Jj  =  cos'1  ^ 


(A.1) 


(A.2)' 


(A.3) 


*  f  Vi 

maximum  relative  elevation  (see  Fig.  A.2):  emx  =  (3j  -  cos-1  — * - 

L  Tj  -  Ar/2  J 

r  COS0„  -  COS(3-  COSP 

maximum  relative  azimuth  (see  Fig.  A.3):  jxmx  =  cos'1  - $ — — ^ - 

sinPj  sin(3 


Figure  A.l  Footprint  center, range  bounds, and  jth  range  cell. 
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2 


Figure  A.2  Maximum  relative  elevation  emx  ( A  =A  r/2). 
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APPENDIX  B 

-4 

Relevant  Information  for  a  Scatterrr  C-  (see  Fig.  B.l). 

4  ■*  * 

B.l  Computation  of  {r  *^r. from  [m,C-,a  } 

4  4  4  4 

(m  ,C  •)  arc  the  vertical  projections  of  (m,Cp 

ground  range:  GR  =  |  C^  -  mg  |  =  J  [(Xj  -  dx)2  +  -  dy)2] 

range  to  C  •  =  r  •  =  |  C  ■  -  m  |  =  J  [GR2  +  h2} 

6  gi  gi  1  gi  ' 

slant  range  to  C-:  rj  =  |  C.  —  m  |  =  J  [GR2  +  (h  -  h|)2] 

*  .* 

range  cell  index:  j  =  INT[(r-  —  r  +  Ar/2)/Ar]  +  j 

-l  f  GR 

elevation  angle:  P-  =  tan  — ^ — 

azimuth  angle:  =  tan_1[(Xj  -  dx)/(yi  -  dy)] 

* 

relative  azimuth  increment:  (J.-  =  a-  -  a 

* 

B.2  Computation  of  {r.a^p^x^}  from  [r  } 

/  2  2 
ground  range:  GR  =  y  (r^  —  h  ) 

slant  range  to  C-:  r-  =  J  [GR2  +  (h  —  h-)2] 

range  cell  index  j  from  r  is  as  above 

* 

azimuth  angle:  a-  =  a  + 

★ 

elevation  angle:  Pj  =  Pj  + 
x— coordinate:  x^  =  GR  sin  cXj  +  dx 
y-coordinate:  y^  =  GR  cos  cr  +  dy 

•4 

B.3  Antenna  Gain  for  C-  at  time  t. 

-4  -4*  -4-4-4 

Antenna  angle  to  C^:  Ang=  Angle(q  -  a(t),  C-  -  a(t)) 

let  wl  =  (ax(t)  ,  ay(t)  ,  az(t)} 

*  *  * 

w2  =  {qx,  qy)qz} 
and  w3  =  {Cix,Ciy)Ciz} 

s  =  sin(Ang)  =  J  (1  -DOTCOS(wl,w2,w3)**2) 


(B.l) 


(B.2) 


(B.3a) 
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a=0  line 
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where  the  function  DOTCOS  is  listed  in  Appendix  E. 

One-way  Antenna  gain  =  DISH(s,2,DOVL)  (B.3b) 

where  the  function  DISH  is  listed  in  Appendix  C,  and  DOVL  is  the  dish  diameter  in 


wavelength  units. 


The  gain  is  computed  at  three  times 

D  ,  D  . 

{t  = - f  ,  t  =  0  ,  t  =  -^} 


(B.3c) 
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APPENDIX  C 

Listings  for  DISH  and  BESSJ1 
These  functions  are  used  to  compute  the  dish  antenna  gain: 


FUNCTION  DISH(SINBET,IDDISH,DDISH) 

C  SINBET=SINE  OF  OFF-BROADSIDE  ANGLE 
EPS=.000001 

IF(ABS(SINBET)XT.EPS)  GO  TO  1 
PI=3. 141592654 
GO  TO  (1,2,3),  IDDISH 

1  DISH=1. 

RETURN 

2  X=SINBET*PI*DDISH 
DISH=ABS(2.*BESSJ1(X)/X) 

RETURN 

3  BET=ASIN(SINBET) 

X=PI*DDISH*BET 
IF(ABS(X).LT.EPS)  X=EPS 
DISH=ABS(SIN(X)/X) 

RETURN 

END 

FUNCTION  BESSJl(X) 

REAL*8  YrPl,P2TP3,P4,P5,Ql,Q2,Q3,Q4,Q5,Rl,R2,R3,R4,R5,R6, 

*  S1,S2,S3,S4,S5,S6 

DATA  R1,R2,R3,R4,R5,R6/72362614232.D0,— 7895059235.D0.242396853.1D0 

* 

*’  -297261 1 .439D0, 1 5704.48260D0  -30. 1 6036606D0/, 

*  S  1,S2,S3,S4,S5,S6/144725228442 £>0,2300535 178.DO, 

*  1 8583304.74D0,99447.43394D0, 376.999 1 397D0, 1  .DO/ 

DATA  P1J*2J>3JP4JP5/1.D0,.183105D— 2,— .3516396496D—4,.2457520174D— 5 

* 

*’  — .240337019D-6/,  Ql,Q2,Q3,Q4,Q5/.04687499995DO,-.2002690873D-3 
* 

*’  .8449199096D— 5,-.88228987D-6,.105787412D-6/ 

EF(ABS(X)£T.8.)THEN 

Y=X**2 

BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))» 

*  /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6))))) 

ELSE 

AX=ABS(X) 

Z=8./AX 

Y=Z**2 

XX=AX— 2.356 19449 1 

BESSJl=SQRT(.636619772/AX)*(COS(XX)*(Pl+Y*(P2+Y*(P3+Y*(P4+Y 

*  *P5))))— Z*SIN(XX)*(Ql+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)s'/)) 

*  *SIGN(1.,X) 

END  IF 
RETURN 

END 
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APPENDIX  D 


Details  of  Taylor  Expansion  of  T- 

This  appendix  shows  the  computation  of  the  Taylor  expansion  given  in  (5.9).  This 
Taylor  expansion  differs  in  three  important  ways  from  that  presented  in  appendix  A  of  [2]  pan 
II:  (1)  it  accounts  for  scatterer’s  height,  (2)  the  notation  is  more  general,  and  (3)  the  derivation 
is  more  elegant. 

The  range  r-  is  expressed  as  a  function  of  four  incremental  variables  {h,t,|i,£}  in  (4.4). 
Our  strategy  is  to  expand  r-,  then  the  round  trip  delay  r.  is  obtained  through  multiplication  by 
2/c.  The  notation  for  the  partial  derivatives  was  explained  in  (4.8),  where  an  upper  bar  means 
evaluation  for  zero  increments.  In  our  application  the  incremental  variables  are  defined  by 

<vl’v2’v3’v4}  =  CD-I) 

2 

It  is  convenient  to  introduce  a  new  variable  y  =  r  - , 
y  =  h2  +  h2  tan2((3+e)  +  v2t2  -  2hvt{tan(P+£)sinpv  cos(a+p.-av)-cos(3v }  (D.2) 


where  (3  stands  for  p  • ,  a  stands  for  a  ,  and  r  stands  for  r.  in  this  appendix. 

The  expression  for  y  is  simplified  by  introducing  two  new  functions,  f  and  g: 


y  =  h2  +  h2!2  +  v2t2  -  2hvt{fg  —  cospy} 

(D.3a) 

f  =  tan(P  +  £) 

(D.3b) 

g  =  cos(a  +  \i  -  ay)  sinpy  . 

(D.3c) 

The  chain  rule  is  used  to  compute  the  partial  derivatives. 

the  partial  derivatives  of  r  with  respect  to  y. 

It  follows  that  we  need  all 

ry  =  0.5y-°-5 

(D.4a) 

ryy  *  -°.25y-1-5 

(D.4b) 

ryyy=0.375y-2'5  . 

Evaluation  for  zero  increments  yields 

(D.4c) 
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f  =  tan($ 

(D.5a) 

g  =  cos(a  -  ay)  sin(3y 
y  =  h2  +  h2  f 2 

(D.5b) 

(D.5c) 

II 

1  U 

(D.5d) 

ry  =  0.5/r 

(D.5e) 

f yy  =  -0-25/f3 

(D.5f) 

iyyy  =  0.375/f5  . 

(D.5g) 

The  expressions  for  the  partial  derivatives  of 

f  and  the  corresponding  reference 

values 

are 

f  -  f  -  ^ 

?  _  1 

(D.6a). 

4  e  cos2(P+e) 

4  cos2]} 

f  ^  2sin((5+e) 

44  66  cos3(J3+e) 

f  2sing 

44  3d 

cos  p 

(D.6b) 

c  f  2  +  4sin2(0+e) 

444  “  cos4(|3+e) 

?  2  +  4sin23 

r444  4f» 

cos  p 

(D.6c) 

Similarly,  the  expressions  for  the  partial  derivatives  of  g  and  the  corresponding  reference 

values  are 

g3  =  gp  =  -sin(a+p-av)sinPy 

g3  =  -sin(a-av)sinpy 

(D.7a) 

g33  =  g^  =  -cos(a+p-ay)sinpv 

g33  =  -cos«x-<xv)sin(3v 

(D.7b) 

g333  =  gppp  =  sin(a+n-av)sinpy 

g333  =  sin(a-<Xv)sinPv  • 

(D.7c) 

The  expression  for  the  partial  derivatives  of  y  and  the  corresponding  reference  values 


are 
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y.  =-2H-2Vjt 
y2  =  2v“t  —  fg  +  2Hv. 

y3  -  -%■ 

y4  =  2ff8-fe«t 

i'll"2 

*12  "  "2vz 
*13  =0 
*14  =  0 

*  22  * 2y2 

*23  *  “f  V 
*24  -  -feg  =  -f4« 

*33  "  "fW 
*34  *  -feV 
*44 y 2{l*21 
*111 
*112  =  0 
*113  =  0 
*114  ~  0 
*122  =  0 
*123  =  0 
*124  =  0 

*133  =° 

*134=° 

*144  =  ° 

*222  m  0 
*223  -  0 


yj  =  — 2h 

y2  =  -f  |  +  2hvz 

y3  =  o 

*4  -  2f  f4 
*11  =2 
*12  *  -2vz 
*13  =° 

*14  =  ° 

*22  -  2*2 
*23  =  _f  ®3 
*24  -  -f4* 

*33  =  0 
*34  =  ° 

*44  =  2F4  +  2f  f44 
*111=° 
yn2  =  0 
yn3  =  ° 

yu4  =  o 
*122  =  0 
*123  =0 
*124  =  0 
*133  =  0 
^134  =  0 
yi44=0 
*222 =  0 
*223  =  0 


(D.8) 
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y224  =  0 

y224  =  0 

y233  =  _fV 

y233  =  ^33 

y234  =  s\i 

y234  =  -F4  % 

y244  = 

y244  =  ~F44® 

y333  =  ~ 

y333  =  0 

y334  ” 

y334  =  0 

y344  ^ee^jx* 

y344  =  0 

y444  -  4fe  ^££+  2fE  ^E£+  2^  ^££E  ^EEE^ 

y444  =  6 F4F44 

Finally  all  the  reference  values  of  the  partial  derivatives  of  range  r  are: 

fl*fy*l  '  f2  =  fy  ^2  '  F3=0  •  f4”fy*4 
'll’ Vll  +  fyy*l 
f  12  *  fy  ^12  +  fyy  *1  ?2 
f13  =0 

'l4  “  'yy  ^1  ^4 

f22  *  fy  5-22  +r‘yy  *2 

r23  =  ry  y23 

?24  =  Fy  y24  +  Fyy  y2  y4 

r33  =  ° 

r34  =  0 

r44  _  ry  y44  +  ryy  y4 
r  1 1 1  “  ^ryy  yl  yll  +  ryyy  yl 
f112  -  'yy^l  S' 12  +  h  i’ll)  +  fyyy  h 

rH3  =° 

r  1 14  ~  ryy  y4  yll  +  ryyy  yl  y4 

r"l22  =  Fyy^2y2  y12  +  yly22>  *  r’yyy  y2  yl 


(D.9) 
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r123  =  ryy  yl  y23 

r"l24  =  ?yy(yl  y24  +  y4  y12>  +  r’yyy  yl  y2  y4 

r133=0 

r134  =  0 

fj44  =  ryy  yi  ^44  +  *yyy  y4  yl 
^222  =  ^ryy  y2  y22  +  ^yyy  y2 


r223  “  2ryy  y2  y23 

?224  =  ?yy(2y2  y24  +  y4  y22>  +  fyyy  y2  y4 
r233  =  ry  y233 


r234  =  ry  y234 


y4y23 


r244  =  ry  y244  +  ryy^2y4  y24  +  y2  y44^  +  ryyy  y4  y2 
r333  =  0 
r334  =  0 


f344  =  0 

* aaa  =  y^  +  3r. 


y 4  y44 


The  Taylor  expansion  for  range  r  is  given  by 

.2  _  t2  _ 

r(l, 2,3.4)  =  r+(r1h,+r2t+r4e)  +  {i’ll  J"  +  ^12^  +  ^14^  +  r22  2  +  r23tJA 


+  r24t£  +  r44  2 


c2  - 
7  }+trm  5 


h2t  -  h2e  -  ht2  - 
+  r  1 12  T  +  r  1 14  ~T  +  r  122  ~T  +  r123 


+  r^^hlli+r.-j^fUE 


2  3  2  2  ^  2  _  -te2 

+  r144  ~T  +  r222  ST  +  r223  2  +  r224  X  +  r233  +  r234  ^  +  r244  ~2~ 


r444  5“' 


(D.10) 


For  the  round  trip  delay,  multiply  (D.10)  by  2/c.  This  is  the  same  as  replacing  ?  by  x,  r  by 
x,  by  x  and  r by  x^.  where  i,  j,  k  take  the  values  { 1,2,3,4}.  The  Taylor  expansion  for 
the  round  trip  delay  is  in  (4.9) 


68 


APPENDIX  E 


Listings  for  FRRDP  and  Associated  Subprograms 
This  Appendix  contains  the  listings  for  FRRDP  and  for  the  associated  subprograms  in 
alphabetical  order.  Excluded  are:  DISH  and  BESSJ1  which  are  in  Appendix  C. 


PROGRAM  FRRDP  .FOR 


C  ASSUME  NBSCA  IN  GROUND  MAP  AND  A  FOOTPRINT  AT  WQ9  ,1/16/89,  1700  hr 
SLARGE 

CHARACTER*240  TITLE 

REAL  W(27),TIMCNT,TIMDWL,VELOMvALFM,BETM,WMPOSC(3),WQ9(5),DELRAN, 
&FFNBIT,WPAR(6)AIFXMT,SELWIN,BETKAI1SELDIS,BITINT 
REAL  SMOMAP(450)  ,TXD(64)  ,TXQ(64) , 

&WADRLO(165),WADRUP(165),WADSMO(90),WAGAIN(3),WANGN9(165), 

&WAGNCL(  128),WBET9(  1 65),  WDUM(3),WDO  VL(3),WEPSMX(  1 65),  WIND(64) 

REAL  WMPOSB(3),WMPOSE(3),WMU3DB(165),WMUMDR(166),WMVELO(3), 

&WRAN9(  1 65),  WTAU  (25)  ,WTIME(64),WTSMO(5),  WW  APEX(  1 28 ,9), 

&XD(  128),XM(128),XQ(  1 28) 

INTEGER  LNPPP(128) 

REAL  DISHDOTCOS,FRQOFM,FREQ 

REAL  CLIGHT JDELFIN,DELFUT,DELTIM,DOVL JEPS,FMU,FREQO,GDISH2, 
&GDISH3,GRC,HDELRA,HTHETA,HTIMDW,RANCOF,RANMN,RANMX,OMEGO,SCAFAC, 
&SINANG,TEMP,THETA,TMPBET,TPI,V1,V2,V3,X»XSCA,YSCA,Y,Z,FZMX 
INTEGER  EDDISH, IFTLE,INDX,INTCOF,IW, IRAN ,JRANQ9JRANMX,JW,K, 
&KADSMO,KDQ,KI,KJ,KFF,KNPPPvKSCALO,KSCAUP,KT,KTIM,KTHSCA, 

&nbitff,nbitfx,nfft,nffth,nffxvnffxh^npppmn,npppmx,nscafp,nxmap, 

&NYMAP 

REAL  SYNMAP(128,55) 

EQUIVALENCE  (W(l),TIMCNT),(W(2),TIMDWL),(W(3),VELOM),(W(4)ALFM), 
&(W(5)3HTM),(W(6),WMPOSC(  1  )),(W(9),WQ9(  1  )),(W(  1 5),DELRAN),(W(  1 6), 
&FFNBIT),(W(17),WPAR(l)),(W(23)^aFXMT),(W(24),SELWIN),(W(25)3ETKA 
&U,(W(26),SELDIS),(W(27)3ITINT) 

C  PARAMETERS  FOR  OUTPUT  MAP  ARE  IN  WPAR=(DELMAP,THRESHXMNMAP,YMNMAP 

C  WQ9*CXQ3cQ3QALFQ9,BETQ93ANQ9)  ********************************* 

C  DIMENSION  OF  SMOMAP  =  NBSCA*  15  (SEE  SUB.  CLUTINFO)  *************** 

C  INPUT  FILES:  CRRDPI.INP  AND  CLUTINFI.INP  ************************* 

C  OUTPUT  FILES:  FRRDPG.GRD  AND  FRRDPP.PRT  ************************** 

FRQ0FMU(V1,FMU,V2,V3)=-V1*(V2*FMU+(V3/2.)*FMU**2) 

0PEN(UNIT=113ILE=,CLUTINFI.INF,STATUS=,0LD') 

OPEN(UNIT=9,FILE=’CRRDPI.INF,STATUS=’OLD') 

OPEN  (UNIT =3  ,FTLE='FRRDPP.  PRT  ,ST  ATU  S=’OLD') 
OPEN(UNIT=7vFILE='FRRDPG.GRD'!STATUS='OLD') 
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C  HALF  3db  BEAMWIDTH=.0215  gives  DO VL=  14.83343 168=DISH  DIAMETER  IN 

C  &LAMBDA  UNITS  (  2-WAY  GAIN  =  -3.00  DB  means  DISH=0.84139514)******* 

DATA  THETA,FREQO,CLIGHT,TPI,EPS,IFILE,WDOVL/.043,94.E9,3.E8, 

&6.283 1 8531, .0001, 3, 17. 15094, 17. 15094, 14.80884/ 

SINANG=SIN  (0.02 1 5) 

GDISH2=DISH(SINANG,2,WDOVL(2)) 

GDISH3=DISH(SINANG,3,WDOVL(3)) 

WRTIE(3,*)’GDISH2,GDISH3,,GDISH2,GDISH3 

TITLE='TIMCNT=1;  TIMDWL=2;  VELOM=3;  ALFM=4;  BETM=5;  WMPOSC=6-8;  WQ 
&9=9-l0  DUMMY=11— 14;DELRAN=15;FFNBIT=16;WPAR1=17=DELM 

&AP- WPAR2=  1 8=THRESH;WPAR3= 1 9=XMN  ;WPAR4=20=YMN/  WPAR5=21=XMX;  WPAR6 
&=22=YMX;  XIFXMT=23,SELWIN=24;  BETKAI=25;SELDIS=26;  BmNT=27' 

CALL  WINPUT(27,W,9  .TITLE) 

CALL  PRTXY(W,  1  ,27,28, IFILE) 

C  READJUST  (XMNMAP.YMNMAP)  TO  BE  INTEGER  MULTIPLE  OF  DELMAP  ******** 

DO  3  K=3  4 

3  WPAR(K)=WPAR(  l)*FLOAT(INT  (WPAR(K)/WPAR(  1 ))) 

EDDISH=tNT(SELDIS+.00001) 

DOVL=WDOVL(IDDISH) 

NXMAP=  1  .+(WPAR(5)— WPAR(3))/WPAR(  1) 

NYMAP=  1  .+<WPAR(6)— WPAR(4)  )/WPAR(  1 ) 

WRITE(*,*),NXMAP=,JDCMAP/  NYMAP=',NYMAP 

OMEGO=TPI*FREQO 

HDELRA=DELRAN/2. 

HTHETA=THETA/2. 

NBITFF=INT(FFNBIT+EPS) 

NBITFX=NB  ITFF+INT  (B ITINT +EPS ) 

NFFT=2**NBITFF 

NFFTH=NFFT/2 

NFFX=2**NBITFX 

NFFXH=NFFX/2 

DELFIN=  1 VTIMD  WL 

INTCOF=2*  *  INT(B  ITINT+EPS ) 

DELFUT =DELFIN/FLO  AT  (INT  COF) 

DELTIM=TIMDWL/NFFT 
DO  2  KTIM=1,NFFT 

2  WTIME(KTEM)=DELTIM*(KTIM-NFFrH-0.5) 

C  COMPUTE  MISSILE  INFORMATION  ************************************** 
WMVELO(l)=VELOM*SIN(BETM)*SIN(ALFM) 

WMVELO(2)= VELOM*S  IN(BETM)  *COS  ( ALFM) 

WMVELO(3)=— VELOM*COS(BETM) 

HTIMDW=TIMDWL*.5 
DO  21  K=l,3 

WMPOSB(K)=WMPOSC(K)+WTIME(  1  )*  WMVELO(K) 

21  WMPOSE(K)=WMPOSC(K)+WTIME(NFFT)*WMVELO(K) 

C  COMPUTE  FOOTPRINT  INFORMATION  AT  WQ9  ***************************** 
GRC=SQRT((WQ9(  1  }-WMPOSC(  1  ))**2+(WQ9(2)-WMPOSC(2))**2) 
WQ9(4)=ATAN2((WQ9(l)-WMPOSC(l)),(WQ9(2)-WMPOSC(2))) 
WQ9(5)=ATAN2(GRC,WMPOSC(3)) 

WQ9(6)=SQRT(GRC**2+WMPOSC(3)**2) 

RANMN=WMPOSC(3)/COS(WQ9(5)-HTHETA) 

RANMX=WMPOSC(3)/COS(WQ9(5)+HTHETA) 

JRANQ9=INT((WQ9(6)-RANMN+HDELRA)/DELRAN) 
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JRANMX=2*JRANQ9 

C  STORE  RANGE, BET,EPSMX,AND  MU-3DB  FOR  FICTITIOUS  SCATTERERS  ******* 

DO  7  JRAN=1,JRANMX 

WRAN9(JRAN)=WQ9(6)+(JRAN— JRANQ9)*DELRAN 
WBET9(JRAN)=ACOS(WMPOSC(3)/WRAN9(JRAN)) 
TMPBET=ACOS(WMPOSC(3)/(WRAN9(JRAN)-HDELRA)) 
WEPSMX(JRAN)=WBET9(JRAN)— TMPBET 
7  WMU3DB(JRAN)=ACOS((COS(HTHETA)-COS(WQ9(5))*COS(WBET9(JRAN)))/ 
&(SIN(WQ9(5))*SIN(WBET9(JRAN)))— 0.00000016) 

C  INITIALIZE  ARRAY  SYNMAP  ****************************************** 

DO  13  JRAN=1  JRANMX 
DO  13  KNPPP=1,NFFX 
13  SYNMAP(KNPPPJRAN)=0. 

C  THIS  IN  PREPARATION  FOR  MAPPING  ********************************** 

C  ENTER  CLUTTER  INFORMATION  FOR  FOOTPRINT  AT  WQ9  AND  MAP  PARAMETERS 
WRITE(3,*)'JRAMNX=' JRANMX 

CALL  CLUTINFO(DELRANjrrHETA,JRANQ9 JRANMX,WQ9,WMPOSC,WPAR,NSCAFP, 
&SMOMAP,WADRLO,WADRUP,WADSMO) 

C  INDEX  KNPPP=1,NFFX  FOR  FREQ=[-NFFXU*DELFUT,(NFFXH— 1  )*DELFUT]  ***** 

C  COMPUTE  ARRAY  LNPPP,THEN  KFF=LNPPP(KNPPP)  ************************ 

DO  9  KNPPP=  1  ,NFFX 
INDX= KNPPP+NFFXH 
IF(INDX.GT.NFFX)  INDX=INDX-NFFX 
9  LNPPP(KNPPP)=INDX 

C  PROCESS  RANGE  CELLS  SEQUENTIALLY  AND  COMPUTE  REFLECTIVITY  MAP  **** 
DO  35  JRAN=1  JRANMX 
KSCALO=WADRLO(JRAN)+EPS 
EF(KSCALO.EQ.O)  GO  TO  35 

C  INITIALIZING  QUADRATURE  ARRAYS=(XDJCQ)  *************************** 

DO  11  K=1,NFFX 
XD(K)=0. 

11  XQ(K)=0. 

WDUM(3)=0. 

C  COMPUTE  TAYLOR  COEFFICIENTS  WTAU=(TAU,TAUt,TAUe,TAUtt,TAUee,TAUtm, 

C  TAUte,TAUttt,TAUeee,TAUtmrn,TAUtcc,TAUtme,TAUttm,TAUtte)  ********** 

CALL  TAYLOR(VELOM,BETMALFM,WQ9(4),WQ9(5),WMPOSC(3),HTHETA, 
&DELRANJRAN.WTAU) 

WRirE(3,*)’ARRAY  WTAU 
CALL  PRTXY(WTAU,  1 ,25, 119,3) 

C  PREPARATION  FOr’mAPPING:  COMPUTE  (NPPPMN,NPPPMX)  ON  3DB  FOOTPRINT 
C  TO  LIMIT  #  CELLS  BEING  MAPPED  TO  KNPPP=NPPPMN,NPPPMX  ************* 
NPPPMN=FRQOFMU(FREQO,WMU3DB(JRAN),WTAU(6),WTAU(10))/DELFUT 
&+NFFXH— 5. 

NPPPMX=FRQOFMU(FREQO,— WMU3DB(JRAN),WTAU(6),WTAU(10))/L'ELFUT 
&+NFFXH+5. 

C  COMPUTE  MU  AT  CENTER  OF  CELLS  FOR  INCREASING  DOPPLER  ************* 

DO  17  K=NPPPMN— l.NPPPMX+l 
FREQ=DELFUT*FLOAT(K-NFFXH-l)/FREQO 
17  WMUMDR(K)=(— WTAU(6)+SQRT(WTAU(6)**2-2.*WTAU(10)*FREQ))/WTAU(10) 

C  COEFFICIENT  FOR  RANGE  NORMALIZATION:  RANCOF  ********************** 
RANCOF=(WRAN9(JRAN)AVRAN9(JRANQ9))**4 
KSCAUP=WADRUP(JRAN)+EPS 
DO  30  KTHSCA=KSCALO,KSCAUP 
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KADSMO=WADSMO(KTHSCA)+EPS 

C  INFORMATION  ON  SCATTERER,WTSMO=(AMPSCArFASSCA,MUSCAJRAN-SCA,ZSCAj 

DO  31  KT=lr5 

31  WTS  MO(KT)=S  MOMAP(KADS  MO— 4+KT) 

C  ANTGAIN(t)=WAGAIN(l)+WAGAIN(2)*t+WAGAIN(3)*tA2  ******************* 

C  WAGAIN  IS  OUTPUT  OF  SCAGAN, INPUTS  ARE  DOVLrHTIMDW,WMPOS=(BEG, 

C  &CENTJBND),WQ9,WTSMO=(AMPJ7AS,MUrJRAN,ZSCA)JRANQ9,DELRAN  ********* 
CALL  SCAGAN (IDDISHJDO  VLJITIMD  W,  WMPOSB ,  WMPOS  C, WMPOSE , W Q9 , WTS  MO, 
&JRANQ9J5ELRAN, WAGAIN) 

CALL  TSIGNAL(OMEGO,NFFT,RANCOF,XIFXMT,WTSMO,WTAU)WTIME,W  AGAIN, TXD, 

&TXQ) 

DO  33  KDQ=1,NFFT 
XD(KDQ)=XD(KDQ)+TXD(KDQ) 

33  XQ(KDQ)=XQ(KDQ)+TXQ(KDQ) 

30  CONTINUE 

CALL  DERAMP(OMEGO,WTAU(4),NFFT,XD,XQ,WTIME) 

CALL  WINDOWS(NFFT,4,BETKAI,WIND) 

DO  37  K=1,NFFT 
XD(K)=XD(K)*WIND(K) 

37  XQ(K)=XQ(K)*WIND(K) 

CALL  FFT(-NBITFX,  1,1,1  ,XD,XQ) 

C  COMPUTE  ARRAY  OF  ANTENNA  GAIN  AT  CENTER  OF  CELLS  ***************** 

DO  39  KNPPP=NPPPMN,NPPPMX 

SINANG=SQRT(1.— (SIN(WQ9(5))*SIN(WBET9(JRAN))*COS(WMUMDR(KNPPP))+ 
&COS(WQ9(5))*COS(WBET9(JRAN)))**2) 

39  WAGNCL(KNPPP)=DISH(SINANG,IDDISHJX)VL)**2 

C  (XD,XQ)  NEED  TO  BE  SCALED  FOR  RANGE  AND  ANTENNA  GAIN  AT  KNPPP  **** 

DO  41  KNPPP=NPPPMN,NPPPMX 
KFF=LNPPP(KNPPP) 

SCAFAC=  1 7(RANCOF*  WAGNCL(KNPPP)) 

XD(KFF)=XD(KFF)*SCAFAC 
XQ(KFF)=XQ(KFF)*SCAFAC 
XM(KFF)=SQRT  (XD(KFF)  *  *2+XQ(KFF)  *  *2) 

C  SYNMAP  STORES  THE  FOOTPRINT  MAP,  X<=>JRAN;  Y<=>DOPPLER  ********* 
IF(XM(KFF).GT.  1 .0)  SYNMAP(KNPPP,JRAN)=XM(KFF) 

WRITE(3,*)7RAN=’,FLOAT(JRAN),'  KNPPP=',FLOAT(KNPPP),’  SYNMAP=’, 
&SYNMAP(KNPPPJRAN) 

41  CONTINUE 
35  CONTINUE 
FZMX=0. 

DO  60  IW=1,NFFX 
DO  60  JW=1  JRANMX 

EF(SYNMAP(IWJW).GT.FZMX)  FZMX=SYNMAP(IWJW) 

60  CONTINUE 

WRITE(7;rDSAA”)') 

WRITE  (7,' (15, 1X,I5)')NFFX  JRANMX 
WRITE(7;  (E 1 2.5, 1XE 1 2.5)’)  1 .  ELOAT(NFFX) 

WRITE(7,’{E  1 2.5, 1XE 1 2.5)’)  1  .ELOAT(  JRANMX) 

WRITE(7,’(E  1 2.5, 1XJE 1 2.5)')0.  EZMX 
DO  69  JW=1  JRANMX 

WRITE(7,62)  (SYNMAP(IW,JW),IW=I,NFFX) 

WRITE(7,’( )’) 

69  CONTINUE 
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62  FORMAT(10(1XF8.3)) 
CL0SE(3) 

CL0SE(9) 

CL0SE(7) 

CLOSE(ll) 

STOP 

END 


SUBROUTINE  CLUTINFOEOR 


SUBROUTINE  CLUTINFO(DELRAN>HTHETA,JRANQ9JRANMX,WQ9,WMPOSC1WPAR, 
&NSCAFP,SMOMAP,WADRLO,WADRUP,WADSMO) 

CHARACTER*240  TITLE 

REAL  WMPOSC(l),WSCA(7),WQ9(l),WPAR(l),WMSCA(3),SMOMAP(l), 
&GRDMAP(150),WADRLO(1),WADRUP(1),WJRAN(200),WADSMO(1),W(34) 

REAL  Wl,COSANT,XSCA,YSCA,HDELRA,DOTCOS,EPS,RGI,XI,YI,FMUI 
INTEGER  NW,KK,KW,NBSCA,KSMO,KSCA,KXJRANON,KSCAT,KINC,KJ 
C  WQ9=(XQ,YQ2Q^LFQ9,BETQ9,RANQ9)  ********************************* 

C  DIMENSION  OF  GRDMAP  IS  NBSCA*4  *********************************** 

C  FIRST  CALL  ENTER  THE  GROUND  MAP  INTO  GRDMAP  ********************** 

C  INPUTS  FOR  GRDMAP:NBSCA,SCAAMP,SCAFAS,SCAX,SCAY  ****************** 
TITLE=' 

&  ENTERING  SCAT-DATA,  W(1)=3*NUMBER  OF  SCAT+1, 

&  { W (3i— 1  ),W(3i),W(3i+ 1 ) } = { Xi,  YiJHfoari }  or  {MUi,RGi,Hbari}  forith 
&SCATTERER' 

READ(11,2)  W1 
2  FORMAT(F10.5) 

NW=INT(ABS(W 1 )+.000 1 ) 

REWIND  11 

CALL  WINPUT(NW,W,  1 1  .TITLE) 

COS  ANT=COS(HTHETA*  1.13) 

KK--4) 

DO  3  KW=2,NW,3 
GRDMAP(KK+ 1  )=  1 . 

GRDMAP(KK+2)=0. 

XI=W(KW) 

YI=W(KW+1) 

FMUI=ATAN2((XI— WMPOSC(l)),(YI— WMPOSC(2)))— WQ9(4) 

RGI=SQRT((XI— WMPOSC(l))**2+(YI— WMPOSC(2))**2+WMPOSC(3)**2) 

EFOVl.GT.O.)  GO  TO  1 
FMUI=W(KW) 

RGI=W(KW+1) 

XI=SQRT(RGI**2-WMPOSC(3)**2)*SIN(FMUI+WQ9(4))+WMPOSai) 
YI=SQRT(RGI**2— WMPOSC(3)**2)*COS(FMUI+WQ9(4))+WMPOSC(2) 

W(KW)=XI 

W(KW+1)=YT 

1  WRITE(*,99)XI,  YI,  FMUI,  RGI,  W(KW+2) 

99  FORMAT(1X,,XI=,E10.5,1X,’YI=,E10.5,1X,’MUI=\E12.5,1X;RGI=,> 

&E10.5,1X,'HI=’E10.5) 

GRDMAP(KK+3)=W(KW) 
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GRDMAP(KK+4)=W(KW+ 1 ) 

GRDMAP(KK+5)=W(KW +2) 

3  KK=KK+5 

NBSCA=(NW— 1)/3 

C  SMOOTHED  OUTPUT  FOR  CURRENT  FOOTPRINT  IS  SMOMAP,  ONE  REAL 
C  &SCATTERER  GIVES  15-COMPONENTS  ************************************ 

C  SMOMAP=(EARLY,SCAFAS,MUJRANON— l,ZSCA,ON,SCAFAS,MU>JRANON,ZSCA,LATE 
C  &,SCAFAS^IU JRANON+ 1  ,ZSCA)  ***************************************** 

KSMO=0 

HDELRA=DELRAN/2. 

C  WSCA=(X-SCA,Y-SCA,Z-SCA,AMP-SCA,FAS-SCA)  ************************* 

DO  5  KSCA=1,NBSCA 
KX=(KSCA— 1)*5+3 
W  SCA(  1)=GRDMAP(KX) 

WSCA(2)=GRDMAP(KX+1) 

WSCA(3)=GRDMAP(KX+2) 

W  SCA(4)=GRDMAP(KX— 2) 

WSCA(5)=GRDMAP(KX— 1) 

IF(DOTCOS(WMPOSC,WQ9,WSCA).LT.COSANT)  GO  TO  5 
EF(W5CA(1).LT.(WPAR(3)— WPAR(l)).OR.WSCA(l).GT.(WPAR(5)+WPAR(l))) 

&GOTO5 

IF(WSCA(2).LT.(WPAR(4)— WPAR(l)).OR.WSCA(2).GT.(WPAR(6)+WPAR(l))) 

&GO  TO  5 

C  REPLACE  (X,Y)  BY  (MU.RANSCA)  ************************************* 
WSCA(6)=WSCA(1) 

WSCA(7)=WSCA(2) 

WSCA(  1  )=AT  AN2(WSCA(6)— WMPOSC(  1  ),WSCA(7)-WMPOSC(2))-WQ9(4) 
WSCA(2)=SQRT((WSCA(6)-WMPOSC(l))**2-KWSCA(7)-WMPOSC(2))**2+ 

&(WSCA(3>— WMPOSC(3))**2) 

C  RANGE  SMOOTHING:  1-SCATTERER  GIVES  WMSCA=(EARLY,ON,LATE)  ********** 
JRANON=((WSCA(2)-WQ9(6)+HDELRA)/DELRAN)+JRANQ9 
WRITE(3,*),JRANON=,JRANON;  JRANQ9=',JRANQ9 

CALL  RANGSMOO(DELRAN,JRANONJRANQ9,WQ9(6),WMPOSC,WSCA,WMSCA) 

KSCAT=0 

DO  7  KINC=2,12,5 

SMOMAP(KSMO+KINQ=WSCA(5) 

SMOMAP(KSMO+KINC+ 1  )=WSCA(  1 ) 

KS  CAT=  KS  CAT+ 1 

SMOMAP(KSMO+KINC— 1)=WMSCA(KSCAT) 

SMOMAP(KSMO+KINC+2)=JRANON— 2+KSCAT 
7  SMOMAP(KSMO+FQNC+3)=WSCA(3) 

KSMO=KSMO+ 1 5 
5  CONTINUE 
EPS  =0.00001 
NSCAFIM) 

DO  10  K=4,KSMO,5 
NSCAFP=NSCAFP+ 1 
WADSMO(NSCAFP)=K 
10  WJRAN(NSCAFP)=SMOMAP(K) 

CALL  SORTRI(WJRAN,WADSMO ,-NSCAFP) 

DO  12  J=1,JRANMX 
WADRLO(J)=0. 

WADRUP(J)=0. 


74 


DO  14  K=NSCAFP,1,— 1 

14  IF(ABS(WJRAN(K)-FLOAT(J)).LE.EPS)  WADRLO(J)=FLOAT (K) 
DO  16  K=lrNSCAFP 

16  IF(ABS(WJRAN(K)-FLOAT(J)).LE.EPS)  WADRUP(J)=FLOAT(K) 
12  CONTINUE 
RETURN 
END 


SUBROUTINE  DERAMP JPOR 


SUBROUTINE  DERAMP(OMEGO,TAU4,NFFT,XD,XQ,WTIME) 
REAL  WTIME(1),XD(1),XQ(1) 

DOUBLE  PRECISION  HOMTTT,TIMtANGROT,TXD,TXQ 

HOMTTT=0.5*OMEGO*TAU4 

DO  5  K=1,NFFT 

TIM=WTIME(K) 

ANGROT =HOMTTT*TIM*  *2 
TXD=XD(K) 

TXQ=XQ(K) 

XD(K)=TXD*DCOS(ANGROT)— TXQ*DSIN(ANGROT) 
XQ(K)=TXD*DSIN(ANGROT)+TXQ*DCOS(ANGROT) 

5  CONTINUE 
RETURN 
END 


FUNCTION  DOTCOSJFOR 


FUNCTION  DOTCOS(Wl,W2,W3) 

C  DOTCOS=COS{ANG(W2— W1  W3— 1)}  ************************************** 

REAL  W1(1),W2(1),W3(1),V1(3),V2(3), VIM, V2M 
V1(1)=W2(1)— Wl(l) 

V 1  (2)=W2(2)— W 1  (2) 

V 1  (3)=W2(3)— W  1(3) 

V1M=SQRT(V1(1)**2+V1(2)**2+V1(3)**2) 

V2(1)=W3(1)— Wl(l) 

V2(2)=W3(2)— W  1(2) 

V2(3)=W3(3)— Wl(3) 

V2M=SQRT(V2(1)**2+V2(2)**2+V2(3)**2) 

DOTCOS=(V  1  ( 1  )*’ V2(  1 )+ V 1  (2)*V2(2)+V1(3)*V2(3))/(V  1 M*  V2M) 
IF(ABS(DOTCOS).GE.  1 .)  DOTCOS=SIGN(0.99999999,DOTCOS) 

RETURN 

END 


SUBROUTINE  FFTFOR 
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SUBROUTINE  FFT(L,KRZIN,KRZOUT,KUNSCR^R,XI) 

REAL  XR(l),XI(l),TPI*EPSFLEPSMO,HOLDR,ANG,TR,TI,DELANG,HOLDI 
INTEGER  M,NJ,I,K  J 1 JG, JUMP, JUMPPO  JUMPMO J 1  MAX, JUMPEX,KMAX,IMAX, 
&NDISTM,KANDS,NDIST,NSPANrNHALF,NHALFP 
TPI=6.283 185307 17958 
M=IABS(L) 

EPSFI=l.E-30 
EPSMO=l.E— 15 
EF(M.EQ.O)  RETURN 
N=2**M 

GO  TO  (9,10),KRZIN 
10  DO  25  J=1,N 
HOLDR=XR(J) 

XR(J)=HOLDR*COS(XI(J)) 

25  XI(J)=HOLDR*SIN(XI(J)) 

9  NHALF=N/2 

NHALFP=NHALF+1 
DELANG=TPI/FLOAT  (N) 

EFflLXT.O)  DELANG=-DELANG 

NSPAN=N 

NDIST=NSPAN/2 

DO  2  1=1  >1 

ANG=0.0 

DO  3  J=1,NDIST 

TR=COS(ANG) 

TT=SIN(ANG) 

DO  5  K=J,N,NSPAN 

KANDS=K+NDIST 

HOLDR=XR(K)-XR(KANDS) 

HOLDI=XI(K)-XI(KANDS) 

XR(K)=XR(K)+XR(KANDS) 

XI(K)=XI(K)+XI(KANDS) 

XR(KANDS)=TR*HOLDR-TI*HOLDI 
XI(KANDS)=TI*HOLDR+TR*HOLDI 
5  CONTINUE 
3  ANG=ANG+DELANG 
NSPAN=NDIST 
NDIST=NSPAN/2 

2  DELANG=DELANG+DELANG 
IF(L.LT.O)  GO  TO  23 
GO  TO  (18,19),KRZOUT 

18  DO  1  J=1,N 

XR(  J)=XR(  J)/FLOAT  (N) 

1  XI(J)=XI(J)/FLOAT(N) 

GO  TO  21 

19  DO  20J=1,N 
HOLDR=XR(J)/FLOAT(N) 

HOLDI=XI(J)/FLOAT(N) 

XR(  J)=SQRT  (HOLD  R*  *2+HOLDI*  *2) 

IF(XR(J).GT.EPSMO)  GO  TO  26 

XI(J)=0.0 

GO  TO  20 
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26  IF(ABS(HOLDR)XT.EPSFI)  HOLDR=SIGN(EPSFI,HOLDR) 
XI(J)=ATAN2(H0LDIrH0LDR) 

20  CONTINUE 

21  CONTINUE 
GO  TO  4 

23  GO  TO  (4,ll),KRZOUT 

11  DO  22  J=1,N 
HOLDR=XR(J) 

XR(J)=SQRT(HOLDR**2+XI(J)**2) 

EF(XR(J).GT.EPSMO)  GO  TO  27 

XI(J)=0.0 

GO  TO  22 

27  IF(ABS(HOLDR)XT.EPSFI)  HOLDR=SIGN(EPSFIJIOLDR) 
XI(J)=ATAN2(XI(J),HOLDR) 

22  CONTINUE 

4  IF(M.EQ.l)  RETURN 

EF(KUNSCR.EQ.O)  RETURN 

JUMPEX=1 

NDIST=N 

IMAX=M/2 

DO  6  I=1,IMAX 

JUMP=JUMPEX 

JUMPPO= JUMP+ 1 

JUMPMO=JUMP— 1 

JUMPEX=JUMPEX+JUMPEX 

NSPAN=NDIST 

NDIST=NDIST/2 

NDISTM=NDIST— JUMP 

DO  8  JG=JUMPPO,N,NSPAN 

J1MAX=JG+NDISTM 

DO  8  J1=JG,J1MAXJUMPEX 

KMAX= J 1  +  JUMPMO 

DO  8  K=J1,KMAX 

KANDS=K+NDISTM 

HOLDR=XR(K) 

HOLDI=XI(K) 

XR(K)=XR(KANDS) 

XI(K)=XI(KANDS) 

XR(KANDS)=HOLDR 

8  XI(KANDS)=HOLDI 

6  CONTINUE 
RETURN 
END 


SUBROUTINE  K3NPUT.FOR 


SUBROUTINE  KINPUT(K) 
INTEGER  K 
WRirEI^.'GX/^'W)') 
READ(*,’(I10)’)K 
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RETURN 

END 


SUBROUTINE  PAYOFXEOR 


SUBROUTINE  PAY0FX(X1,X2,X3,Y1,Y2,Y3,A1A2A3) 

REAL  X12M^C12P^C13M,X13P 
X12M=X1— X2 
X12P=X1+X2 
X13M=X1— X3 
X13P=X1+X3 

A3=((Y1— Y2)*X13M— (Yl— Y3)*X12M)/(X12M*X12P*X13M— X13M*X13P*X12M) 

A2=((Y  1— Y2)— A3  *X  1 2M*X  1 2P)/X  1 2M 

A  1=Y  1-A2*X  1  — A3*X  1  *X  1 

RETURN 

END 


SUBROUTINE  PRTXY  J=OR 


SUBROUTINE  PRTXY (XY ,KST ART,KEND , ISEQ,IFILE) 

REAL  XY(1) 

INTEGER L 

1  FORMAT  (6(  1X,E  1 2. 6)) 

2  FORMAT(lX,  13HPRTXY  KSTART=,I8,18X,5HKEND=,I8,8X,5fflSEQ=J4) 

3  F0RMAT(1X35H***********************************) 
WRITE(IFILE,3) 

WRITE(IFILE,2)KSTART>KEND,ISEQ 

WRITE(IFILE,1)(XY(L)L^KSTART,KEND) 

RETURN 

END 


SUBROUTINE  RANGSMOO  JOR 


SUBROUTINE  RANGSMOO(DELRANJRANONJRANQ9,RANGE9,WMPOSC,WSCA,WMSCA) 
REAL  WMPOSC(l),WSCA(  1  ),WMSCA(  1  )E^C,RANON  ,RANER,RANLT,GRANON, 
&GRANSCJDELON  DHLER  JDELLT.SUM 
INTEGER  K 
F(X)=1.-ABS(X)/1.5 

RANON=RANGE9+(JRANON-JRANQ9)*DELRAN 

RANER=RANON-DELRAN 

RANLT =RANON+DELRAN 

GRANON=SQRT(RANON**2-WMPOSC(3)**2) 

GRANSC=SQRT(WSCA(2)**2-WMFOSC(3)**2) 

DELON=(GRANON-GRANSQ/DELRAN 
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DELER=(DELON-DELRAN)/DELRAN 
DELLT=(DELON+DELRAN)/DELRAN 
WMSC  A(  1  )=F(DELER) 
WMSCA(2)=F(DEL0N) 
WMSCA(3)=F(DELLT) 

SUM=WMSCA(  1  )+WMSCA(2)+WMSCA(3) 
DO  5  K=1  3 

5  WMSCA(k)=WSCA(4)*WMSCA(K)/SUM 
RETURN 
END 


SUBROUTINE  SCAGANJFOR 


SUBROUTINE  SCAGAN(IDDISHJX)VLOTMDW)WMPOSB,WMPOSC,WMPOSE,WQ9, 
&WTSMO,  JRANQ9  JDELRAN,W  AGAIN) 

C  ANTGAIN(t)=WAGAIN(l)+WAGAIN(2)*t+WAGAIN(3)*tA2  ******************* 

C  WAGAIN  IS  OUTPUT  OF  SCANTGAN  .INPUTS  ARE  ©DISHJX)VL,HTIMDW,WMPOS 
C  &=(BEG,CENT^ND),WQ9,WTSMO=(AMP,FASNtU,JRAN,ZSCA),JRANQ9,DELRAN  *** 
C  WGAIN=S CAT— ANTENNA  GAIN  AT  (BEG.CENT.END)  OF  DWELL  *************** 
REAL  DOVL,HTIMDW,WMPOSB(l),WMPOSC(l),WMPOSE(l),WQ9(l),WTSMO(l), 
&DELRAN,WAGAIN(l)1RANSCA,GRANSC,WSCA(3),WGAIN(3),SINANGJX)TCOSJ)ISH 
INTEGER  IDDISHJRANQ9 

RANSCA=WQ9(6)+(WTSMO(4)-FLOAT(JRANQ9))*DELRAN 

GRANSC=SQRT(RANSCA**2-<WMPOSC(3)-WTSMO(5))**2) 

WSCA(l)=WMPOSC(l)+GRANSC*SIN(WQ9(4)+WTSMO(3)) 

WSCA(2)=WMPOSC(2)+GRANSC*COS(WQ9(4)+WTSMO(3)) 

WSCA(3)=WTSMO(5) 

SINANG=SQRT(  1  .-DOTCOS(WMPOSB,WQ9,WSCA)**2) 
WGAIN(l)=DISH(SINANG,IDDISILDOVL)**2 
SINANG=SQRT(  L-DOTCOS(WMPOSC,WQ9,WSCA)**2) 
WGAIN(2)=DISH(SINANG,IDDISH,DOVL)**2 
SINANG=SQRT(  1  .-DOTCOS(WMPOSE,WQ9,WSCA)**2) 
WGAIN(3)=DISH(SINANG,IDDISH,DOVL)**2 

CALL  PAYOFX(-HTIMDW,0. OTMDW,WGAIN(  1),WGAIN(2),WGAIN(3),WAGAIN(  1 ) 
&,WAGAIN(2),WAGAIN(3)) 

RETURN 

END 


SUBROUTINE  SORTRLFOR 


SUBROUTINE  SORTRI  (VALXOC.NS) 
REAL  VAL(1),F,XS 
INTEGER  LOC(l)^,K,MJ,I,n 
N=IABS(NS) 

IF(NS.LT.O)  GO  TO  7 
DO  6  K=1,N 
6  LOC(K)=K 
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7  M=N 

1  CONTINUE 
M=M/2 

IF(M.EQ.O)  RETURN 

K=N-M 

J=1 

2  CONTINUE 
I=J 

3  CONTINUE 
n=I+M 

IF(VAL(DXT.VAL(II))  GO  TO  4 
F=VAL(I) 

XS=LOC(I) 

VAL(I)=VAL(ID 
LOC(I)=LOC(ID 
VAL(II)=F 
LOC(ID=XS 
1=1— M 

IFa.GE.l)  GO  TO  3 

4  CONTINUE 
J=J+1 

IF(J.GT.K)  GO  TO  I 
GO  TO  2 
RETURN 
END 


SUBROUTINE  TAYLORPOR 


SUBROUTINE  TAYL0R(VEL0M,BHTM^JLF\tALFQ3HTQ>HIGHT,HTHETArDRJCELL, 
&WTAU) 

C  HBAR<=>1,  T<=>2,  MU<=>3,  EPS<=>4,  TAU<=>X  ******************* 

C  WTAU=X,X2,X4>X22,X44,X23,X24,X222,X444,X233  /  X244.X234.X223 JC224, 

C  &X1,X11,X12,X14,X111,X112  /  X114,X122,X123,X124,X144  ************** 

REAL  WTAU(  1  ).C0F.RQ,RMN30,VEL0M3ETNLALFM,ALFQ,BETQJHIGHT, 
&HTHETAJ)R3ETA0343443444,G,G3,Y,Y1,Y2,Y4,Y11,Y12,Y22,Y23,Y24, 

&  Y  44,  Y233,  Y234.Y244,  Y444.R.RY  ,R  YY.R  YYY,  VZ.F 
INTEGER  JQJCELL 
VZ*— VELOM*COS(BETM) 

COF=273.0E8 

RQ=H3GHT/COS(BETQ) 

RMN=fflGHT/COS(BETQ-HTHETA) 

JQ=ENT((RQ-RMN+DR/2.)/DR) 

ro=rq+(jcell-;q)*dr 

BETAO=ACOS(fflGHT/RO) 

F=fflGHT*TAN(BETAO) 

F4=HIGHT/COS(BETAO)**2 

F44=2.*fflGHT*SIN(BETAO)/COS(BETAO)**3 

F444=2.*HIGHT*(l.+2.*SIN(BETAO)**2)/COS(BETAO)**4 

G=2.*VEL0M*SIN(BETM)*C0S(ALFQ-ALFM) 

G3»-2.*VELOM*SIN(BETM)*SIN(ALFQ-ALFM) 
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Y=F**2+fflGHT**2 

R=Y**5 

RY=.5/R 

RYY=-.25/R**3 

RYYY=0.375/R**5 

Yl=-2.*fflGHT 

Y2=-E*G+2.*HIGHT*VZ 

Y4=2.*F*F4 

Yll=2. 

Y12=-2.*VZ 
Y22=2.  *  VELOM**2 
Y23=-F*G3 
Y24=— F4*G 

Y44=2.*F4**2+2.*F*F44 

Y233=F*G 

Y234=^?4*G3 

Y244=—F44*G 

Y444=6.*F4*F44+2.*F*F444 

WTAU  ( 1  )=COF*R 

WTAU(2)=COF*RY*Y2 

WTAU(3)=COF*RY*Y4 

WTAU(4)=COF*(RY*Y22+RYY*Y2**2) 

WTAU(5)=COF*(RY*Y44+RYY*Y4**2) 

WTAU(6)=COF*RY*Y23 

WTAU(7)=COF*(RY*Y24+RYY*Y2*Y4) 

WTAU(8)=COF*(3.*RYY*Y2*Y22+RYYY*Y2**3) 

WTAU  (9)=C0F*(R  Y*  Y444+3.  *R  YY*  Y  4*  Y  44+RYYY*  Y  4**3) 

WTAU  ( 10)=COF*R  Y*  Y233 

WTAU(ll)=COF*(RY*Y244+RYY*(2.*Y4*Y24+Y2*Y44)+RYYY*Y4**2*Y2) 

WTAU(12)=COF*(RY*Y234+RYY*Y4*Y23) 

WTAU  ( 1 3)=COF*2.  *RYY*  Y2*  Y23 

WTAU(14)=COF*(RYY*(2.*Y2*Y24+Y4*Y22)+RYYY*Y2**2*Y4) 

WTAU(  15)=C0F*R  Y*  Y1 
WTAU(16)=C0F*(RY*Y1 1+RYY*Y1**2) 

WTAU(17)=C0F*(RY*  Y12+RYY*  Y 1  *  Y2) 

WTAU  ( 1 8)=C0F*R  YY*  Y 1  *  Y4 
WTAU(19)=C0F*(3.*RYY*Y1*Y11+RYYY*Y1**3) 

WTAU(20)=COF*(R  YY*(2.  *  Y 1  *  Y 1 2+ Y2*  Y 1 1  )+R  YYY*  Y 1  **2*  Y2) 
WTAU(21)=C0F*(RYY*Y4*Y11+RYYY*Y1**2*Y4) 
WTAU(22)=COF*(RYY*(2.*  Y2*  Y 1 2+Y 1  *  Y22)+R  YYY*  Y2**2*  Y 1 ) 
WTAU(23)=C0F*RYY*Y1  *Y23 

WTAU(24)=COF*(RYY*(Y  1  *  Y24+ Y4*  Y 1 2)+R  YYY*  Y 1  *  Y2*  Y4) 
WTAU(25)=C0F*(RYY*Y1*Y44+RYYY*Y4**2*Y1) 

RETURN 

END 


SUBROUTINE  TSIGNALJOR 


SUBROUTINE  TSIGNAL(OMEGO,NFFT,RANCOFrXIFXMT,WTSMO,WTAU,WTIME, 
&WAGAIN,TXD,TXQ) 
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REAL  WTSMO(  1  ),WTAU(  1  ),WTIME(  1  ),WAGAIN  ( 1  ),TXD(  1  ),TXQ(  1 ), TAMP  AMP 
DOUBLE  PRECISION  PSI,TPI,TPSI,COFO,COFI,COF2,TIM 
INTEGER  KF 
TPI=6.283 185307 

COFO*-OMEGO*(WTAU(  1)+WTAU(  15)*WTSMO(5))+XIFXMT 
COF1=-OMEGO*(WTAU(6)*WTSMO(3)+.5*WTAU(10)*WTSMO(3)**2+WTAU(17)* 
&WTSMCK5)) 

COF2==-0.5*OMEGO*(WTAU(4)+WTAU(13)*WTSMO(3)) 

TAMP=WTSMO(l)*RANCOF 

DO  5  KF=1,NFFT 

TIM=WTIME(KF) 

AMP=TAMP*(WAGAIN(1)+WAGAIN(2)*TIM+WAGAIN(3)*TIM*TTM) 

TPSI=COFO+COF1  *TIM+COF2*TTM**2 

PSI=DMOD(TPSI,TPI) 

TXD(KF)=AMP*DCOS(PSR 
TXQ(KF)= AMP*DSIN  (PSD 
5  CONTINUE 
RETURN 
END 


SUBROUTINE  WCHANGEPOR 


SUBROUTINE  WCBANGE(W) 

REAL  W(l) 

INTEGER K 

1  PRINT*, ENTER  INDEX  OF  W  TO  CHANGE;  0  MEANS  NO  CHANGE  ' 
CALL  KINPUT(K) 

EF(KLEQ.O)  RETURN 
CALL  XINPUT(W(K)) 

GOTO  1 
END 


SUBROUTINE  WINPUTEOR 


SUBROUTINE  WINPUT(NW,W,IDFTLE, TITLE) 
CHARACTER*240  TITLE 
CHARACTER*55  WCHAR(51) 

REAL  W(l) 

INTEGER  KWJC 
DO  7  KW=1,NW 

7  READ(IDFILE,11)  W(KW),WCHAR(KW) 

11  FORMAT(F16.7A55) 

PRINT*,’ ARRAY  W  BEFORE  UPDATE' 
WRITE(*,’(1XA240)’)  TITLE 
WRITE  (*,21  (W(K),K=1,NW) 

2  FORMAT(4(2XFl  6.7)) 

CALL  WCHANGE  (W) 
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REWIND  IDFILE 
DO  8  KW=1,NW 

8  WRITE(EDFILE,'(F16.7,A55)')  W(KW),WCHAR(KW) 

CLOSE(IDFILE) 

RETURN 

END 


SUBROUTINE  XINPUTJOR 


SUBROUTINE  XENPUT(X) 
REAL  X 

wRrrE(*,'(ix;,x=?",\)') 

READ(*/(F16.7)')X 

RETURN 

END 
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APPENDIX  F 

Listings  for  POLONFRA  and  Associated  Subprograms 
This  Appendix  contains  the  listings  for  POLONFRA  and  its  associated  subprograms, 
except  for  SORTRI  which  is  in  Appendix  E. 

SUBROUTINE  POLONFRA-FOR 


SUBROUTINE  POLONFRA(NBA, WAX, WAY,WPAR,UOJJO,IMXJMX,TAREA,WW AREA, 
&IOKRO.IFAIL) 

C  (WAX,WAY)=ABSOLUTE  IN  METERS,  (XMN,XMX,nOJJO,ZDIST)= ABSOLUTE  IN 
C  DELMAP  UNITS, (A,WX,WY,WG)=RELATTVE  IN  DELMAP  UNITS  *************** 
REAL  WAX(  1),WA  Y  ( 1),WAXT  (5),  WA  YT(5),WPAR(  1  )VA(5,2),ZDIST(54), 

&WWAREA(0: 15,0: 15)X>HORI(5)JD  VERT(5)T>SIDE(5),WG(2),WXRD(54),WYRD 
&(54),EPS,EPSH,OLDDST 

INTEGER  LOC2L 1  (54)  JLOCCL 1  (54)  JLOCRLC(54)  JRECW  (54),IDROTW  (16), 

&NBCINR(0: 15,0: 15)  JHVAW(54),INTCOR(0:16,0: 16)JCONTZ,IDREC 
COMMON  /FIRST/KONTZ,INTCOR,WXRD,WYRDJHVAW,EPS^EPSH  /SECOND/OLDDST 
&,WGT>SIDE,DHORI,DVERT 

REAL  AREA,AREAPO,FIG,FJG,SUMX,SUMY,TAREA,TXMN,TYMN,XMED,XMN,XMX, 
&YMED,YMN,YMX 

INTEGER  I,IC,  IMXM,INEXT,IOKRO  J  J C JMXMJL-KONT,  K,KONT C,K  1 ,  K2T  C, 

&LCPJLR.L 1  ,L  1  A,L 1 B  JL2,NRJN4 
EXTERNAL  AREPOL 
DO  1  K=1,NBA 

WAXT(K)=(W AX(K)— WPAR(3))/WPAR(  1 ) 

1  WAYT(K)=(W  A  Y(K)-WPAR(4))/WPAR(  1 ) 

C  IF  CONVEX  &  COUNTERCLOCKWISE  ROTATION  ,IOKRO=0,  THEN  SKIP  CONROT  * 
EF(IOKRO.NE.O)  CALL  CONROT(NBA,WAXT,WAYT,IFAIL) 

EF(IFAIL.EQ.l)  GO  TO  55 
EPS=  0001 
EPSH=EPS/2. 

XMN=WAXT(1) 

YMN=WAYT(1) 

DO  5  K2=2,NBA 
TXMN=WAXT(K2) 

TYMN =W  A YT  (K2) 

IF(TXMN1T  JCMN)  XMN=TXMN 
IF(TYMNUT.YMN)  YMN=TYMN 
5  CONTINUE 

IIO=INT(XMN-0.0001) 

JJO=INT(YMN-0.0001) 

C  TRANSFER  THE  NBA-APEXES  INTO  A(3,2)  IN  RELATIVE  DELMAP  UNITS  ***** 

DO  3  K1=1,NBA 

A(K1 , 1  )=WAXT(K1  )-FLOAT(HO) 

3  A(K1 ,2)=WA  YT(K1 )— FLOAT(JJO) 
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non 


TAREA=AREPOL(NBA,WAXT,WAYT) 

IF(TAREAJLT.EPS)  RETURN 

C  COMPUTE  UPPER  MAP  INDEXES  FOR  FRAME  :  (IMXJMX)  ****************** 
XMX=A(1,1) 

YMX=A(U) 

DO  7  K2=2,NBA 

IF(A(K2, 1).GT.XMX)  XMX=A(K2,1) 

EF(A(K2,2).GT.YMX)  YMX=A(K2,2) 

7  CONTINUE 
IMX=ENT(XMX+.99999) 

JMX=INT(YMX+.99999) 

IMXM=IMX— 1 
JMXM=JMX— 1 

C  COMPUTE  THE  CENTER  OF  GRAVITY  WG  ********************************* 
SUMX=0. 

SUMY=0. 

DO  8  K2=I,NBA 
SUMX=SUMX+A(K2, 1 ) 

8  SUMY=SUMY+A(K2V2) 

WG(l)=SUMX/FLOAT(NB  A) 

WG(2)=SUMY/FLOAT(NBA) 

C  INITIALIZATION  OF  INTCOR  ***************************************** 

DO  9  K1=0,IMX 
DO  9  K2=0  JMX 

9  INTCOR(K1,K2)=0 

C  INITIALIZE  WWAREA  ************************************************ 

DO  33  I=0,IMXM 
DO  33  J=0,JMXM 
33  WWAREA(U)=0. 

C  FIND  INTERSECTS  AND  INTERIOR  CORNERS  FOR  EACH  SIDE  OF  THE  POLYGON 
C  SEQUENTIALLY  (A1  TO  A2,  A2  TO  A3,  A3  TO  Al,...  A(NBA-l)  TO  A(NBA)) 

KONTZ=0 
OLDDST=0.0 
DO  10  K1=1,NBA 
K2=K1+1 

EF(K2.EQ.(NBA+1))  K2=l 

DSIDE(K  1  )=SQRT  (( A(K  1 , 1 )— A(K2, 1))**2+(A(K1 ,2)—A(K22))**2) 
DHORI(Kl)=ABS(A(Kl,l)-A(K2,l)) 

D  VERT  (K 1  )= ABS  ( A(K  1 ,2)— A(K2,2)) 

C  COMPUTING  WXRD(LI),WYRD(L1),IHVAW(L1),ZDIST(L1)  FOR  Ll=l,KONTZ  ***** 
CALL  UNEAB(IMX  JMX,K1  A(K1 , 1  ),A(K1 ,2)^A(K2, l)w\(K2,2)ZX>IST) 

10  CONTINUE 
IF(KONTZ.GT.NBA)  GO  TO  4 

C  KONTZ=NBA  MEANS  THE  POLYGON  IS  WITHIN  ONE  RECTANGLE  *************** 
T=INT(WG(1)) 

J=ENT(WG(2)) 

WWAREA(U)=TAREA 

RETURN 

4  LKONT=KONTZ 

CALL  SORTRI(ZDISTLOC2L  1  ,LKOND 

AFTER  SORTRI  ZDIST(L2)  INCREASES  WITH  L2  WHERE  L2=l,KONTZ 
AND  L1=L0C2L1(L2)  ************************************************ 

FOR  AN  INTERIOR  CORNER  INTCOR(U)=NBA,  MAKE  IT  1  ****************** 
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DO  25  IC=0,IMX 
DO  25  JC=OJMX 

25  INTCOR(ICJC)=INTCOR(ICJC)/NBA 

C  ROTATE  CC  ON  CONTOUR  AND  MARK  CORNER  POINTS  ********************* 
DO  13  L2=UCONTZ 
L1=L0C2L1(L2) 

IC=INT(WXRD(Ll)+.49999) 

JC=INT  (WYRD(L 1 )+.49999) 

IF(ABS(WXRD(L1)—IC)+ABS(WYRD (LI)— JQ.GT.EPS)  GO  TO  13 
WXRD(L  1  )=FLO  AT  (IC) 

WYRD(L  1  )=FLOAT(JQ 
INTCOR(ICJC)=l 
13  CONTINUE 

C  IF  4  INTERIOR  CORNERS  NBCINR(IJ)=4,THEN  AREA  =  1.  *********** 

DO  23  1=0  JMX 
DO  23  J=0  JMX 

NBCINR(U)=INTCOR(IJ)4-INTCOR(I+ 1  J)+INTCOR(I  J+ 1)+ 

&INTCOR(I+U+l) 

IF(NBCINR(Irr).EQ.4)  WWAREA(I,J)=1. 

23  CONTINUE 

C  CONTOUR  IS  COMPRESSED  FROM  KONTZ  TO  KONTC  BY  COPYING  ONLY  ID  OF  ** 
C  POINT  IFF  CC  DIST  TO  PRIOR  POINT.LT.EPSR  A  COMPRESSED  APEX  FALLS 
C  ON  THE  RECTANGLE  GRID.  OUTPUT  IS  LOCCLl;Ll=LOCCLl(LC)  ************ 
CALL  COMPRESZ(LOC2LlJCONTQLOCCLl) 

C  ROTATE  CC  ON  COMPRESSED  CONTOUR  POINT  ID  IS  LC=1  JCONTC************ 

C  SEGMENT(A(LC),B(LC+1))  IS  BOUND  BY  RECTANGLE  IDREC:IRECW(LC)=IDREC 
C  SKIP  IDREC  IF(A3)  IS  ON  THE  RECTANGULAR  GRID  ******************** 

NR=0 

DO  15  LC=1, KONTC 
LCP=LC+1 

IF  (LCP.GT.KONTC)  LCP=LCP-KONTC 
L 1  A=LOCCL  1  (LQ 
L 1  B=LOCCL  1  (LCP) 

XMED=0.5*(WXRD(L1A)+WXRD(L  1 B)) 

YMED=0.5*(WYRD(L1A)+WYRD(L1B)) 

FIG=FLOAT  (INT  (XMED+0.49999)) 

FJG=FLOAT(INT(YMED+0. 49999)) 

IF(ABS(XMED-nG).LT.EPSRORABS(YMED-FJG).LT.EPSH)  GO  TO  15 

NR=NR+1 

I=INT(XMED) 

J=INT(YMED) 

IDREC=J*IMX+I 

IF(NBCINR(U).EQ.4)  IDREC=IDREC+10(X) 

IRECW  (NR)=IDREC 
LOCRLC(NR)=LC 
15  CONTINUE 

C  SORT  IRECW  FOR  INCREASING  IDREC  .BACKTRACK  ARRAY  IS  LOC3LC  ******* 
CALL  SORTHORECWa-OCRLC ,-NR) 

C  AFTER  SORTII  ERECW(LR)  WHERE  LR=1,NR  AND  LC=LOCRLC(LC)  ******** 

N4=0 

IRECW  (NR+ 1 )=99999 
DO  50  LR=1,NR 

C  IF  RECTANGLE  IDREC  BOUNDS  TWO  SEGMENTS(A  1 B 1  A2B2),EDREC=IDRECW(LR) 
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C  =IRECW(LR+1)VLC(A  1  )=L0CRLC(LR),LC(A2)=L0CRLC(LR+ 1 ).  LCs  FOR  SAME  * 
C  IDREC  ARE  STORED  IN  IDR0TW={LC(A1)J-C(B1)=LC(A1)+1,LC(A2),LC(B2)= 

C  LC(A2)+1 }.  SKIP  IF  IDREC. GE. 1000  (RECTANGLE  ALREADY  MAPPED)  ****** 
IDREC=IRECW  (LR) 

IF(IDREC.GE.  1000)  GO  TO  50 

N4=N4+1 

LC=LOCRLC(LR) 

IDROTW  (N4)=LC 

N4=N4+1 

LCP=LC+1 

IF(LCP.GT.KONTQ  LCP=LCP-KONTC 

EDROTW(N4)=LCP 

INEXT=IRECW(LR+1) 

IF(INEXT.EQ.IDREC)  GO  TO  50 
I=MOD(IDRECJMX) 

J=INT  (IDREC/IMX) 

C  AFTER  SORTI  IDROTW(L5)  INCREASES  WITH  LR  (INCREASING  LCs  )******** 
CALL  SORTI(IDROTW,N4) 

CALL  RECTAREa(UN4,KONTCO-OCCL1, IDROTW  AREA) 

WWAREA(IJ)-AREA 

N4=0 

50  CONTINUE 
RETURN 

55  WRirE(6,*)'POLYGON  NOT  CONVEX,  OR  3  POINTS  ON  ONE  LINE  !!!!!!!!!!’ 
RETURN 
END 


FUNCTION  AREPOLFOR 


FUNCTION  AREPOL(NBX,Y) 

REAL  AREA ATRI AF^CR,XS  A 1 ,X2,X3,X(  1),YF,YR,  YS,  Y 1 , Y2, Y3, Y(  1 ) 
INTEGER  K,NBM 

ATR1(X1,Y1A2,Y2A3,Y3)=0.5*ABS(X1*Y2+X2*Y3+X3*Y1-Y1*X2~Y2*X3- 

&Y3*X1) 

IF(NB.LE.2)  GO  TO  3 
AREA=0. 

XR=X(1) 

YR=Y(1) 

NBM=NB— 1 
DO  5  K=2JNBM 
XF=X(K) 

YF=Y(K) 

XS=X(K+1) 

YS=Y(K+1) 

AREA=AREA+ATRI(XR,YRAF,YF,XS,YS) 

5  CONTINUE 
AREPOL=AREA 
RETURN 

3  PRINT*, ERROR  IN  ARE  POL:  NBA<3' 

RETURN 
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SUBROUTINE  COMPRESZJOR 


SUBROUTINE  COMPRESZ(LOC2L  1  ,KONTCXOCCL  1 ) 

INTEGER  LOC2L  1(1  ),LOCCL  1(1  ),IHV  A  W  (54)  JNTCOR(0: 1 6,0: 1 6),KONTZ,L  1 , 
&L1MJL2J-2M 

REAL  WXRD(54),WYRD(54)EPS^EPSHJDIFDIS 

COMMON  /FIRST/KONTZJNTCOR,WXRD,WYRDJHVAW,EPSEPSH 

KONTC=0 

L2M=KONTZ 

L 1  M=LOC2L  1  (L2M) 

DO  3  L2=lrKONTZ 
L1=L0C2L1(L2) 

DIFDIS= AB  S(WXRD(L  1 )— WXRD(L  1  M))+ ABS  (WYRD(L  1 )— WYRD(L  1 M)) 

IF(DEFDIS.GT.EPSH)  GO  TO  5 

EF(IHVAW(L1)XE.2)  LOCCL 1  (KONTC)=LOC2L  1  (L2) 

GO  TO  7 

5  KONTC=KONTC+ 1 

LOCCL  1  (KONT  C)=LOC2L  1  (L2) 

7  L2M=L2 

L 1  M=LOC2L  1  (L2M) 

3  CONTINUE 
RETURN 
END 


SUBROUTINE  CONROTJOR 


SUBROUTINE  CONROT(NBA,WAX,WAYJFAIL) 

REAL  WAX(  1  ),WAY  ( 1),WAXT(5),  WA  YT(5),T 

INTEGER  K,KONCC,KONCW ,MA,MB,MC 

IFAIL=0 

KONCW=0 

KONCC=0 

DO  1  K=1,NBA 

MA=MOD((K— 1  ),NB  A)+ 1 

MB=MOD(K,NB  A)+ 1 

MC=MOD((K+ 1 )  ,NB  A) + 1 

T=(WAX(MQ— WAX(MB))*(W  A  Y(MA)-WA  Y  (MB))— (WAX(MA)-WAX(MB))*(WA  Y  (MC) 
&— WAY(MB)) 

DF(TJEQ.O.)  GO  TO  3 
IF(T.LT.O.)  KONCW =KONCW + 1 
IF(T.GT.O.)  KONCC=KONCC+ 1 
1  CONTINUE 

IF(KONCC.NE.O-AND.KONCW.NE.O)  GO  TO  3 
IF(KONCC.GT.O)  RETURN 
DO  4  K=1,NBA 
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WAXT(K)=WAX(K) 

4  WAYT(K)=WAY(K) 

DO  5  K=2,NBA 
WAX{K)=WAXT(NBA+2— K) 

5  WAY(K)=WAYT(NBA+2-K) 
RETURN 

3  DFAIL=1 
RETURN 
END 


SUBROUTINE  INSERT JFOR 


SUBROUTINE  INSERT(K,KP,XR,YR,NBINS,AX,AY) 

REAL  AX(1)^Y(1)^R(1),YR(1),DISTAJ>ISTJ)ISTB 
INTEGER  IDISTA,IDISTB,KD,KDT 

C  COMPUTE  EPl=LDROTW(K),IP2=LDROTW(KP)  ON  TRIANGLE  CONTOUR  COUNT 
C  INTERIOR  CORNERS  COUNTERCLOCKWISE  ,#=  NB INS, STORE  RELATIVE  X,Y 
C  INTO  AX,AY.  COORDINATES  OF  IP1,IP2  ARE  IN  XR,YR  ***************** 
NBINS=0 

CALL  XYDIS(XR(K),YR(K)JDISTA,  1) 

CALL  XYDIS  (XR(KP), YR(KP)  JDISTB ,  1 ) 

EF(DISTB.LT.DISTA)  DISTB=DISTB+4. 

IDISTA=INT  (DISTA+ 1 .0000 1 ) 

ID  IS TB=INT(DISTB— 0.00001) 

IF(IDISTB.LT.IDISTA)  RETURN 
DO  5  KD=IDISTA,IDISTB 
KDT=MOD(KD,4) 

DIST=FLOAT(KDT) 

NBINS=NBINS+1 

CALL  XYDIS(AX(NBINS),AY(NBINS),DIST,2) 

5  CONTINUE 
RETURN 
END 


SUBROUTINE  LJNEAB.FOR 


SUBROUTINE  LINEAB(IMX,JMX,K1,AX,AY3X3Y,2X>IST) 

INTEGER  IHVAW(54),INTCOR(0: 16,0: 16)30NTZ,IM JBJCORJ^Q JNO, 

&JJAJBJCOR 

REALDHORI(5)3VERT(5)3>SIDE(5)^DIST(l),WG(2),WXRD(54),WYRD(54), 

&EPS3PSH,0LDDST,0NUN3T0A,PXT0AX3YT0AYREFSIG,TSIGN,TX,TY,UI,UJ, 

&VADI,VADJ,X,Y 

COMMON  /FIRST/KONTZ,INTCOR,WXRD,WYRD,IHVAW,EPS,EPSH  /SECOND/OLDDST 
&,WG,DSIDE3>H0RI,DVERT 

0NLIN(X,YAXAY3X3Y)=(AX-X)*(BY-Y)-(AY-Y)*3X-X) 

C  TESTING  THE  CENTER  OF  GRAVITY  WRT  THE  LINE  A-B  ******************* 
REFSIG=ONLIN(WG(l),WG(2)AXAY3X3Y) 
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C  ENTERING  INFO  ON  APEX  A  ****************************************** 
KONTZ=KONTZ+ 1 
ZDIST(KONTZ)=OLDDST 
WXRD(KONTZ)=AX 
WYRD(KONTZ)=AY 
IHVAW(KONTZ)=7 

C  COMPUTE  UNIT  VECTORS  UJ  ON  Y,AND  UI  ON  X  ************************* 
UJ=(BY— AY)/DSIDE(K1) 

UI=(BX— AX)/DSIDE(K1) 

VADI=0.5+SIGN(0.5,UI) 

VADJ=0.5+SIGN(0.5,UJ) 

IA=INT  ( AX+VADI) 

IB=ENT(BX+1.-VADI) 

JA=INT  (AY +VADJ) 

JB=INT(BY+1.— VADJ) 

DF(ABS(UI).LT.l.E-6)  GO  TO  6 

C  ENTERING  INFO  ON  VERTICAL  INTERSECTS  ***************************** 
INCI=I 

EF(UI.LT.O.)  INCI=-1 
DO  3  I=LA,IBJNCI 
PXTOAX=ABS(I— AX) 

PTOA=DSIDE(Kl)*PXTOAX/DHORI(Kl) 

KONTZ=KONTZ+ 1 

ZDIST  (KONTZ)=PTO  A+OLDDST 

WXRD(KONTZ)=FLOAT(I) 

WYRD(KONTZ)=  AY+UJ*PTOA 
EHVA  W  (KONTZ)=2 
3  CONTINUE 

6  IF(ABS(UJ).LT.  1.E— 6)  GO  TO  10 

C  ENTERING  INFO  ON  HORIZONTAL  INTERSECTS  *************************** 
INCJ=1 

EF(UJ.LT.O)  INa=-l 
DO  5  J=JAJB,INCJ 
PYTOAY=ABS(J— AY) 

PTOA=DSIDE(Kl)*PYTOAY/DVERT(Kl) 

KONTZ=KONTZ+ 1 
ZDIST  (KONTZ)=PTOA+OLDDST 
WXRD(KONTZ)=AX+UI*PTOA 
WYRD(KONTZ)=  FLOAT(J) 

5  IHVAW(KONTZ)=I 
10  CONTINUE 

C  INCREMENT  INTCOR  BY  1  IF  CORNER  (ICORJCOR)  INTERIOR  TO  LINE  AB 
DO  15  ICOR=0,IMX 
DO  15  JCOR=OJMX 
TX=FLOAT(ICOR) 

TY=FLOAT(JCOR) 

TSIGN=0NLIN(TX,TYAXAY,BX3Y) 

IF((TSIGN*REFSIG).GT.O.)  INTCOR(ICOR,JCOR)=INTCOR(ICORJCOR)+ 1 
15  CONTINUE 

OLDDST=OLDDST+DSEDE(K  1 ) 

RETURN 

END 
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SUBROUTINE  RECTAREAJFOR 


SUB  ROUTINE  RECTAREA(I,J,N4,KONTC  J-OCCL 1  ,IDROTW  AREA) 

C  COMPUTE  OVERLAP  AREA  OF  TRIANGLE  &  RECTANGLE  AT  (IJ)  ************ 
INTEGER  EDROTW(  1  ),LOCCL  1  ( l),LDROTW(  1 6)JHVAW(54),INTCOR(0: 1 6,0: 1 6) 

&  JJIVAW(16)  J-DROTX(  14),K4ON,K0N,K0NX,K0NT2  JCPJ-CX 1  ^BINS,NBINSC 
REAL  XR(16),YR(16),XRX(14),YRX(14),AX(4),AY(4),WXRD(54),WYRD(54), 

&AREAPO  JBPS  JEPSH 
EXTERNAL  A  RE  POL 

COMMON  /FIRST/KONTZ,INTCOR,WXRD,WYRD,IHVAW,EPS,EPSH 
C  COMPRESS  EDROTW  INTO  LDROTW  FOR  REPEATED  APEX  (  FROM  N4  TO  KON  )** 
KON=0 
DO  3  K=1,N4 
KP=K+1 

IF(KP.GT.N4)  KP=KP— N4 

IF  (IDROTW (KP).EQ.BDROTW (K))  GO  TO  3 

KON=KON+l 

LDROTW(KON)=IDROTW(K) 

3  CONTINUE 

C  COMPUTE  LOCAL  XY-COORDINATES  IN  (XR,YR),IHVAW  IDENTITY  IN  LHVAW-  ** 
DO  5  K=l,KON 
LC=LDROTW(K) 

Ll-LOCCLl(LQ 

XR(K)= WXRD(L  1 )— FLO  AT  (I) 

YR(K)=WYRD(Ll)-FLOAT(D 
5  LHVAW(K)=IHV  A  W(L  1 ) 

C  EXPAND  ARRAYS  (LDROTW ,XR, YR)  INTO  (LDROTX,XRX,YRX)  FROM  KON  TO 
C  KONX  BY  INSERTING  CORNERS  ON  THE  CONTOUR  ************************* 
KONX=0 
DO  10  K=1JC0N 
KP=K+1 

IF(KP.GT.KON)  KP=KP— KON 
KONX=KONX+ 1 
LDROTX(KONX)=LDROTW  (K) 

XRX(KONX)=XR(K) 

YRX(KONX)=YR(K) 

C  IF  NO  GAP  IN  IDROT  THEN  SKIP  INSERT  ****************************** 
IF((LDROTW(KP)-LDROTW(K)).EQ.1.0R.(LDROTW(KP)-LDROTW(K))JEQ. 
&1-KONTQ  GO  TO  10 
CALL  INSERT(K,KP,XR,YRNBINSVAXVAY) 

DF(NBINS.EQ.O)  GO  TO  10 
DO  15  KIN=1,NBINS 
KONX=KONX+ 1 
LD  ROTX  (KONX)= KIN + 1 00 
XRX(KONX)=AX(KIN) 

YRX(KONX= A  Y  (KIN) 

CONTINUE 

AREA=AREPOL(KONX,XRX,YRX) 

RETURN 
END 
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SUBROUTINE  SORTLFOR 


SUBROUTINE  SORTI  (IVAL^N) 
INTEGER  IVAL(l),M,KJ,I,n,F 
M=N 

1  CONTINUE 
M=M/2 

IFCMEQ.O)  RETURN 

K=N-M 

J=1 

2  CONTINUE 
I=J 

3  CONTINUE 
n=I+M 

DF(IVAL(I).LT.IVAL(II))  GO  TO  4 
F=IVAL(I) 

IYAL(I)=rVAL(II) 

IVAL(II)=F 
1=1— M 

IFa.GE.l)  GO  TO  3 

4  CONTINUE 
J=J+1 

IF(J.GT.K)  GO  TO  1 
GO  TO  2 
RETURN 
END 


SUBROUTINE  SORTILFOR 


SUBROUTINE  SORTH  (VAUXLOCNS) 
INTEGER  VAL(  1  ),XLOC(  1  )fN,K,MJJfII,F,XS 
N=IABS(NS) 

IF(NS.LT.O)  GO  TO  7 
DO  6  K=I»N 

6  XLOC(K)=K 

7  M=N 

1  CONTINUE 
M=M/2 

IF(M-EQ.O)  RETURN 

K=N-M 

J=1 

2  CONTINUE 
I=J 

3  CONTINUE 
n=I+M 

IF(VAL(I).LT.VAL(II))  GO  TO  4 
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F=VAL(D 

XS=XLOC(I) 

VAL(D=VAL(ID 

XLOC(I)=XLOC(II) 

VAL(ID=F 

xloc(ii)=xs 

1=1— M 

EF(I.GE.l)  GO  TO  3 
4  CONTINUE 
J=J+1 

IF(J.GT.K)  GO  TO  1 
GO  TO  2 
RETURN 
END 


SUBROUTINE  XYDIS  .FOR 


SUBROUTINE  XYDIS(X,Y,DIS,KXYDIS) 

INTEGER  ID  IS 

C  KXYDIS=1  FOR  X,Y  TO  DIS,  KXYDIS=2  FOR  DIS  TO  X,Y  ***************** 
GO  TO  (1,2),KXYDIS 

1  IF(X.GE.Y)  DIS=X+Y 
EF(Y.GT.X)  DIS=4.— X— Y 
RETURN 

2  ID  IS =INT(DIS+ 1.00001) 

IF(DIS.LT.O.)  DIS=DIS+4. 

DF(DIS.GT.4.)  DIS=DIS— 4. 

GO  TO  (4,5,6,7),EDIS 

4  X=DIS 
Y=0. 

RETURN 

5  X=l. 

Y=DIS— 1. 

RETURN 

6  X=3.-DIS 
Y=l. 

RETURN 

7  X=0. 

Y=4.— DIS 
RETURN 
END 
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APPENDIX  G 

Listings  for  CRRDP  and  Associated  Subprograms 
This  Appendix  contains  the  listings  for  CRRDP  and  for  the  associated 
subprograms  in  alphabetical  order.  Excluded  are:  the  subprograms  which  are  already 
listed  in  Appendices  C,  D,  E,  and  F. 

PROGRAM  CRRDP .FOR 


C  ASSUME  NBSCA  IN  GROUNDMAP  AND  A  FOOTPRINT  AT  WQ9  ,1/16/89,  1200  hr 
$ LARGE 

CHARACTER*240  'ITT'  E 

REAL  W(27),TIMCNT,TIMDWL,VELOM^\LFM,BETM,WMPOSC(3),WQ9(6)JDELRAN, 
<SJFNBIT,WPAR(6)jaFXMT,SELWINJ3ETKAI,SELDIS,BrriNT 
REAL  SMOMAP(450),TXD(64),TXQ(64), 

&WADRLO(165),WADRUP(165),WADSMO(90),WAGAIN(3),WANGN9(165), 

&WAGNCL(  128),  WBET9(  1 65),  WDUM(3),  WDO  VL(3),  WEPSMX(  165),  WIND(64) 

REAL  WMPOSB(3),WMPOSE(3),WMU3DB(  165),WMUMDR(  166),WMVELO(3), 

&WRAN9(  1 65),WTAU  (25),  WTIME(64),  WTSMO(5),  WW  APEX(  1 28 ,9), 

&XD(  128),XM(128),XQ(  128) 

INTEGER  LNPPP(128) 

RE/*L  DISHJXDTCOSJRQOFM.FREQ 

REAL  CLIGHT,DELFIN,DELFUT,DELTIM,DOVL£PSJ7MU,FREQO,GDISH2, 
&GDISH3,GRC,HDELRA,HTHETA,HTTMDW,RANCOF,RANMN,RANMX,OMEGO,SCAFAC, 
&SINANG,TEMP,THETA,TMPBET,TPI,V1,V2,V3,X.XSCA,YSCA,Y,Z,FZMX 
INTEGER  IDDISHJFDLEJNDX,INTCOF,IWrJRANJRANQ9rIRANMX,JW,K, 
&KADSMO,KDQ,KI,KJ,KFF,KNPPP,KSCALO,KSCAUP,KT,KTIM,KTHSCA, 
&NBITFFJSrarTFXJSlFFT^FFTHJ^FFX^fFFXILNPPPMN^PPPMX,NSCAFP,NXMAP, 

&NYMAP 

REAL  SYNMAP(0:99,0:99),SUMARA(0:99,0:99) 

COMMON  /MAP/S YNMAP.SUMARA 

EQUIVALENCE  (W(l),TIMCNT),(W(2),TIMDWL),(W(3),VELOM),(W(4)ALFM), 

&(W  (5)3ETM),(W  (6),  WMPOSC(  1  )),(W  (9),  WQ9(  1  )),(W(  1 5),DELRAN),(W  (16), 

&FFNBIT),(W(17),WPAR(l)),(W(23)^aFXMT),(W(24),SELWIN),(W(25)3HTKA 

&D,(W(26)3ELDIS),(W(27)3mNT) 

C  PARAMETERS  FOR  OUTPUT  MAP  ARE  IN  WPAR=(DELMAP,THRESH,XMNMAP,YMNMAP 

C  WQ9«(XQ3cQ3Q>ALFQ9,BETQ93ANQ9)  ********************************* 

C  DIMENSION  OF  SMOMAP  =  NBSCA*  15  (SEE  SUB.  CLUTINFO)  *************** 

C  INPUT  FILES:  CRRDPI.INP  AND  CLUTINn.INP  ************************* 

C  OUTPUT  FILES:  CRRDPG.GRD  AND  CRRDFP3RT  ************************** 
FRQOFMU(V  1 3MU.V2,  V3)*-V  1  *(V2*FMU+(V3/2.)*FMU**2) 

OPEN(UNTT=l  l,FILE='CLUTINFI.INP1STATUS=OLD') 

0PEN(UNIT=93ILE=CRRDPI.INF,STATUS=0LD’) 

OPEN(UNIT=*3.FILE-CkRDPP.PRT,STATUS='OLD') 
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OPEN(UNIT=7,FILE='CRRDPG.GRD',STATUS='OLD') 

C  HALF  3db  BEAMWIDTH=.0215  gives  DOVL=  1 4.83343 1 68=DISH  DIAMETER  IN 

C  &LAMBDA  UNITS  (  2-WAY  GAIN  =  -3.00  DB  means  DISH=0.84139514)******* 
DATATHETA,FRE(^,CLIGHT,TPI,EPS,IFILE,WDOVU043,94.E9,3.E8, 

&6.283 1 853 1  ,.0001 ,3, 17. 15094, 17. 15094, 14.80884/ 

SINANG=SIN(0.0215) 

GDISH2=DISH(SINANG,2,WDOVL(2)) 

GDISH3=DISH(SINANG,3,WDOVL(3)) 

WRITE(3,*)'GDISH2,GDISH3',GDISH2,GDISH3 

TITLE='TIMCNT=  1 ;  TIMDWL=2;  VELOM=3;  ALFM=4;  BETM=5;  WMPOSC=6-8;  WQ 
&9=9— 10  DUMMY=11— 14;DELRAN=15;FFNBIT=16;WPAR1=17=DELM 

&AP;WPAR2=18=THRESH;WPAR3=19=XMN  ;WPAR4=20=YMN/  WPAR5=2l=XMX;  WPAR6 
&=22=YMX;  XIFXMT=23,SELWIN=24;  BETKAI=25;SELDIS=26;  BITINT=27' 

CALL  WINPUT(27,W,9,TITLE) 

CALL  PRTXY(W,  1^7,28  JFILE) 

C  READJUST  (XMNMAP.YMNMAP)  TO  BE  INTEGER  MULTIPLE  OF  DELMAP  ******** 

DO  3  K=3,4 

3  WPAR(K)=WPAR(  1  )*FLOAT(INT  (WPAR(K)/WPAR(  1 ))) 

EDDISH=INT(SELDIS+.0000 1 ) 

DOVL=WDOVL(IDDISH) 

NXMAP=1.+(WPAR(5)-WPAR(3))/WPAR(1) 

NYMAP=1.+(WPAR(6)-WPAR(4))/WPAR(1) 

WRrrE(*,*),NXMAP=',NXMAP,'  NYMAP=',NYMAP 

OMEGO=TPI*FREQO 

HDELRA=DELRAN/2. 

HTHETA=THETA/2. 

NBITFF=INT(FFNBIT+EPS) 

NBITFX=NBITFF+INT(BmNT+EPS) 

NFFT=2**NBITFF 

NFFTH=NFFT/2 

NFFX=2**NBITFX 

NFFXH=NFFX/2 

DELFIN=  1  ,/TIMDWL 

INTCOF=2**INT(BmNT+EPS) 

DELFUT=DELFIN/FLOAT(INTCOF) 

DELTIM=TIMDWL/NFFT 
DO  2  KTIM=U^FFT 

2  WTIME(KTIM)=DELTIM*(KTIM-NFFTH-0.5) 

C  COMPUTE  MISSILE  INFORMATION  ************************************** 
WMVELO(l)=VELOM*SIN(BETM)*SIN(ALFM) 

WMVELO(2)= VELOM*SIN  (BETM)  *COS  ( ALFM) 

WMVELO(3)=— VELOM*COS(BETM) 

HTIMDW=TIMDWL*.5 
DO  21  K=l,3 

WMPOSB(K)= WMPOSC(K)+WTIME(  1 )  *  WMYELO(K) 

21  WMPOSE(K)=WMPOSC(K)+WTIME(NFFT)*WMVELO(K) 

C  COMPUTE  FOOTPRINT  INFORMATION  AT  WQ9  ***************************** 
GRC=SQRT((WQ9(  1  )-WMPOSC(  1  ))**2+(WQ9(2)-WMPOSC(2))**2) 
WQ9(4)=ATAN2((WQ9(l)-WMPOSC(l)),(WQ9(2)-WMPOSC(2))) 
WQ9(5)=ATAN2(GRC,WMPOSC(3)) 

WQ9(6)=SQRT(GRC**2+WMPOSC(3)**2) 

RANMN=WMPOSC(3)/COS(WQ9(5)— HTHETA) 


95 


RANMX=WMPOSC(3)/COS(WQ9(5)+HTHETA) 

JRANQ9=INT((WQ9(6)-RA  NMN+HDELRA)/DELRAN) 

JRANMX=2*JRANQ9 

C  STORE  RANGE,BET,EPSMX,AND  MU-3DB  FOR  nCTmOUS  SCATTERERS  ******* 

DO  7  JRAN=1,JRANMX 

WRAN9(JRAN)=WQ9(6)+(JRAN-JRANQ9)*DELRAN 
WBET9(JRAN)=ACOS(WMPOSC(3)/WRAN9(JRAN)) 
TMPBET=ACOS(WMPOSC(3)/(WRAN9(JRAN)— HDELRA)) 
WEPSMX(JRAN)=WBET9(JRAN)— TMPBET 
7  WMU  3DB  (JRAN>= ACOS((COS  (HTHETA)— COS  (WQ9(5))  *COS  (WBET9(  JRAN)))/ 
&(SIN(WQ9(5))*SIN(WBET9(JRAN)  ))-0.000000 1 6) 

C  INITIALIZE  ARRAY  SYNMAP  ****************************************** 

DO  13  IW=0,NXMAP— 1 
DO  13  JW=0,NYMAP— 1 
S  YNMAP(IW  JW)=0. 

13  SUMARA(IW,JW)=0. 

C  THIS  IN  PREPARATION  FOR  MAPPING  ********************************** 

C  ENTER  CLUTTER  INFORMATION  FOR  FOOTPRINT  AT  WQ9  AND  MAP  PARAMETERS 
WRrrE(3,*)'JRAMNX=,JRANMX 

CALL  CLUTINFO(DELRAN,HTHETA,JRANQ9JRANMX,WQ9,WMPOSC,WPAR,NSCAFP, 
&SMOMAP,WADRLO,WADRUP,WADSMO) 

C  INDEX  KNPPP=  1  ,NFFX  FOR  FREQ=[-NFFXH*DELFUT,(NFFXH-1  )*DELFUT]  ***** 

C  COMPUTE  ARRAY  LNPPP.THEN  KFF=LNPPP(KNPPP)  ************************ 

DO  9  KNPPP=  1  ,NFFX 
INDX=KNPPP+NFFXH 
DF(INDX.GT.NFFX)  INDX=INDX-NFFX 
9  LNPPP(KNPPP)=INDX 

C  PROCESS  AND  MAP  RANGE  CELLS  SEQUENTIALLY  ************************* 

DO  35  JRAN=1  JRANMX 
KSCALO=WADRLO(JRAN)+EPS 
EF(KSCALO.EQ.O)  GO  TO  35 

C  INITIALIZING  QUADRATURE  ARRAYS=(XD,XQ)  *************************** 

DO  27  K-LNFFX 
XD(K)=0. 

27  XQ(K)=0. 

WDUM(3)=0. 

C  COMPUTE  TAYLOR  COEFFICIENTS  WTAU=(TAU,TAUt,TAUe,TAUtt,TAUee,TAUtm, 

C  TAUtc,TAUttt,TAUece,TAUtmm,TAUtee,TAUtme,TAUttm,TAUtte)  ********** 

CALL  TAYLOR(VELOM,BETM,ALFM,WQ9(4),WQ9(5),WMPOSC(3),HTHETA, 
&DELRANJRAN,WTAU) 

WRITEO,*)’ ARRAY  WTALT 
CALL  PRTXY(WTAU,1, 25,1 19,3) 

C  PREPARATION  FOR  MAPPING:  COMPUTE  (NPPPMN,NPPPMX)  ON  3DB  FOOTPRINT 
C  TO  LIMIT  #  CELLS  BEING  MAPPED  TO  KNPPP=NPPPMN,NPPPMX  ************* 
NPPPMN=FRQOFMU(FREQO,WMU3DB(JRAN),WTAU(6),WTAU(10))/DELFUT 
&+NFFXH— 5. 

NPPPMX=FRQOFMU(FREQO,-WMU3DB(JRAN),WTAU(6),WTAU(10))/DELFUT 

&+NFFXH+5. 

C  COMPUTE  MU  AT  CENTER  OF  CELLS  FOR  INCREASING  DOPPLER  ************* 

DO  17  K=NPPPMN— 1,NPPPMX+1 
FREQ=DELFUT*FLOAT(K-NFFXH-l  )/FREQO 
17  WMUMDR(K)=(— WTAU(6)+SQRT(WTAU(6)**2— 2.*WTAU(10)*FREQ))/WTAU(10) 

C  COEFFICIENT  FOR  RANGE  NORMALIZATION:  RANCOF  ********************** 
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RANCOF=(WRAN9(JRAN)AVRAN9(JRANQ9))**4 
KSCAUP=WADRUP(JRAN)+EPS 
DO  30  KTHSCA=KSCALO,KSCAUP 
KADSMO=WADSMO(KTHSCA)+EPS 

C  INFORMATION  ON  SCATTERER,WTSMO=(AMPSCAJFASSCA,MUSCAJRAN-SCA,ZSCA) 
DO  31  KT=1  5 

31  WTSMO(KT)=SMOMAP(KADSMO-4+KT) 

C  ANTG AIN (t)=WAG AIN(  1  )+W AGAIN (2) *t+W AG AIN(3)*tA2  ******************* 

C  WAGAIN  IS  OUTPUT  OF  SCAGAN.INPUTS  ARE  DOVLJHTIMDW,WMPOS=(BEG, 

C  «fcCENTJEND),WQ9,WTSMO=(AMPJ:AS,MUJRAN,ZSCA)JRANQ9,DELRAN  ********* 
CALLSCAGAN(IDDISHJXDVLOTMDW,WMPOSB,WMPOSC,WMPOSE,WQ9,WTSMO, 
&JRANQ9,DELRAN, WAGAIN) 

CALL  TSIGNAL(OMEGO,NFFT,RANCOF^GFXMT,WTSMO,WTAU,WTIME,  WAGAIN,  TXD, 
&TXQ) 

DO  33  KDQ=1,NFFT 
XD(KDQ)=XD(KDQ)+TXD(KDQ) 

33  XQ(KDQ)=XQ(KDQ)+TXQ(KDQ) 

30  CONTINUE 

CALL  DERAMP(OMEGO,WTAU(4),NFFT .XD.XQ.WTIME) 

CALL  WINDOWS(64,4,BETKAI,WIND) 

DO  37  K=l,64 
XD(K)=XD(K)*WIND(K) 

37  XQ(K)=XQ(K)  * WIND(K) 

CALL  FFT(-NBITFX,U11VXD,XQ) 

C  COMPUTE  ARRAY  OF  ANTENNA  GAIN  AT  CENTER  OF  CELLS  ***************** 

DO  39  KNPPP=NPPPMN,NPPPMX 

SINANG=SQRT(l.-(SIN(WQ9(5))*SIN(WBET9(JRAN))*COS(WMUMDR(KNPPP))+ 

&COS(WQ9(5))*COS(WBET9(JRAN)))**2) 

39  WAGNCL(KNPPP)=DISH(SINANG,IDDISHJDOVL)**2 

C  (XD,XQ)  NEED  TO  BE  SCALED  FOR  RANGE  AND  ANTENNA  GAIN  AT  KNPPP  **** 

DO  41  KNPPP=NPPPMN,NPPPMX 
KFF=LNPPP(KNPPP) 

SCAFAC=  1 7(RANCOF*WAGNCL(KNPPP)) 

XD(KFF)=XD(KFF)*SCAFAC 
XQ(KFF)=XQ(KFF)*SCAFAC 
41  XM(KFF)=SQRT(XD(KFF)**2+XQ(KFF)**2) 

C  APEXTRAP  COMPUTES  WWAPEX  AT  CENTERS  OF  CEL  1  ,S  ********************* 

C  WWAPEX=(KNPPP;X1,Y1,X2,Y2,X3,Y3,X4,Y4,AREA)  ********************** 

CALL  APEXTRAP(NFFX,WBET9(JRAN),DELRAN,WQ9(4),WEPSMX(JRAN),NPPPMN, 
&NPPPMX,WTAU,WMPOSC,WMUMDR,  WWAPEX) 

C  FFT  CELLS  ARE  MAPPED  EXCEPT  (BELOW  THRESH.or.OUT  FOOTPRT)  KNPPP  ** 

C  GOES : (NPPPMN .NPPPMX)  ie:(MUMX,MUMN);FOR  FFT  USE  KFF=LNPPP(KNPPP)  * 

CALL  MAPRGMCNXMAPJ^YMAP, NPPPMN, NPPPMXrXM,WPARXNPPP, WWAPEX) 

35  CONTINUE 
FZMX=0. 

DO  60  IW=0,NXMAP— 1 
DO  60  JW=0,NYMAP-1 
CF(SUMARA(IW,JW).LT.EPS)  GO  TO  60 
S  YNMAP(TW  JW)=S  YNMAP(IW,JW)/SUMARA(IW  JW) 

IF(SYNMAP(IWJW).GT.FZMX)  FZMX=SYNMAP(IWJW) 

60  CONTINUE 

WRrrE(7,’C'DSAA")') 

WRITE(7,'(I2,1X,I2)')NXMAP,NYMAP 
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WRrrE(7,'(I2,lX,n)')0,NXMAP-l 

WRITE(7,,(I2,1X,I2),)0,NYMAP-1 

WRITE(7,'(I2,1X,I2)’)0,INT(FZMX+EPS) 

DO  69  JW=0,NYMAP-1 

WRTTE(7,62)  (SYNMAP(IW,JW),IW=0,NXMAP-1) 
WRITE(7,'( )’) 

69  CONTINUE 
62  FORMAT(10(1XF8.3)) 

CLOSE(3) 

CLOSE(9) 

CLOSE(7) 

CLOSE(ll) 

STOP 

END 


SUBROUTINE  APEXTRAPJFOR 


SUBROUTINE  APEXTRAP(NPFTUT,BET9T)ELRAN>ALFQ,EPSMX,NPPPMN,NPPPMX, 
&WTAU,WMPOSC,WMUMDR,WWAPEX) 

C  PROGRAM  GENERATED  ON  1/16/1989,  1200  hr  ************************** 

REAL  WTAU(  1 ), WMPOSC(  1  ),WMUMDR(  1  ),WWAPEX(NFFTUT,9),  W(2,2,2),  VARMU, 
&RANMDR,RANMNR,RANMXR,TEMP0,TEMP1 
INTEGER  KTRAP,K1  jC2 

C  TRAPEZOIDS  TO  MAP:KTRAP=(NPPPMN,NPPPMX)  FOR  EACH  KTRAP  COMPUT  MU 
C  AT  (EPS  MN  ,EPS  MX)  FOR  UPPER  BOUNDARY  (UPPER=NEW  LINE,LOWER=OLD  LINE) 
C  FOR  INITIALIZATION  COMPUTE  NEW  LINE  AT  KTRAP=NPPPMN-1  ************ 
VARMU=(WTAU(7)*EPSMX)/WTAU(6) 

RANMDR=WMPOSC(3)/COS(BET9) 

RANMNR=RANMDR-DELRAN/2. 

RANMXR=RANMNR+DELRAN 

TEMP0=— VARMU+(WMUMDR(NPPPMN)+WMUMDR(NPPPMN— 1))*0.5 
TEMP1=TEMP0+2.*VARMU 

W(1,U)=SQRT(RANMNR**2-WMPOSC(3)**2)*SIN(TEMPO+ALFQ)+WMPOSC(1) 

W(2,U)=SQRT(RANMNR**2-WMPOSC(3)**2)*COS(TEMPO+ALFQ)+WMPOSC(2) 

W(l,2^)=SQRT(RANMXR**2-WMPOSC(3)**2)*SIN(TEMPl+ALFQ)+WMPOSC(l) 

W(2,2,2)=SQRT(RANMXR*1*2-WMPOSC(3)**2)*COS(TEMP1+ALFQ)+WMPOSC(2) 

DO  15  KTRAP=NPPPMN,NPPPMX 
DO  10  K1=U 
DO  10  K2=U 

10  W(K1  ,K2, 1  )=W(K  1  ,K2,2) 

TEMP0=s— V  ARMU  -KWMUMDR(  JCTRAP+ 1  )+WMUMDR(KTRAP))*0.5 
TEMP  1  =TEMP0+2.  *  VARMU 

W(1,U)=SQRT(RANMNR*’2-WMPOSC(3)**2)*SIN(TEMPO+ALFQ)+WMPOSC(1) 
W(2,U)=SQRT(RANMNR**2-WMPOSC(3)**2)*COS(TEMPO+ALFQ)+WMPOSC(2) 
W(l,2^)=SQRT(RANMXR**2-WMPOSC(3)**2)*SIN(TEMPl+ALFQ)+WMPOSC(l) 
W(2,2,2)=SQRT(RANMXR**2— WMPOSC(3)**2)*COS(TEMPl+ALFQ)+WMPOSC(2) 
WWAPEX(KTRAP,  1  )= W(  1,1,1) 

WWAPEX(KTRAP,2)=W(2, 1,1) 

WWAPEX(KTRAP,3)=W(  1,2,1) 

WWAPEX(KTRAP,4)=W(2,2,  1 ) 
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WWAPEX(KTRAP,5)=W(  1 ,2,2) 
WWAPEX(KTRAP,6)=W(2,2,2) 
WWAPEX(KTRAP,7)=W(  1,1,2) 
WWAPEX(KTRAP,8)=W(2, 1 ,2) 
15  CONTINUE 
RETURN 
END 


SUBROUTINE  MAPAZMJFOR 


SUBROUTINE  MAPAZM(NXMAP,NYMAP,XMVAL,WAX,WAY,WPAR) 

REAL  WAX(  1  ),WA Y(  1),WPAR(6),WWAREA(0: 1 5,0: 1 5) 

REAL  SYNMAP(0:99,0:99),SUMARA(0:99,0:99) 

REAL  AREPOL 

INTEGER  IJIOJFAIL,IMX,IW,J,JJOJMX,JW 
COMMON  /MAP/SYNMAP,SUMARA 

CALL  POLONFRA(4,WAX,WAY,WPAR,IIOJJO,IMXJMXvAREPOL,WWARjEA,0, 
&IFAIL) 

C  3-SYSTEMS  OF  COORDINATES:  WPAR  define  ABS-location  of  window  ;**** 

C  (IW,JW)=WINDOW— coordinates  ;  (IJ)=FRAME-coordinates  ************* 

DO  43  I=0JMX— 1 
DO  43  J=0,JMX— 1 

DF(WW AREA(I  J).LT.0.002)  GO  TO  43 

iw=no+i 

JW=JJO+J 

IF(IW.LT.O.OR.JW.LT.O)  GO  TO  43 
EF(IW.GT.NXMAP.OR  JW.GT.NYMAP)  GO  TO  43 
SYNMAP(IWJW)=SYNMAP(IW,JW)+WWAREA(U)*XMVAL 
SUMARA(IW,JW)=SUMARA(IWJW)+WWAREA(I,J) 

43  CONTINUE 
RETURN 
END 


SUBROUTINE  MAPRGMPOR 


SUBROUTINE  MAPRGM(NXMAP,NYMAP,NPPPMN,NPPPMX^CM,WPAR,LNPPP,WWAPEX) 
REAL  WPAR(6),WWAREA(0: 15,0: 15), WAX(4),WAY(4),WWAPEX(  128,9), 

&XM(1) 

INTEGER  KA4CAPX,KFF,KNPPPVLNPPP(1) 

DO  45  KNPPP=NPPPMN,NPPPMX 
KFF=LNPPP(KNPPP) 

DF(XM(KFF).LT.WPAR(2))  GO  TO  45 
KA=0 

DO  40  KAPX=1,7V2 
KA=KA+1 

WAX(KA)=WWAPEX(KNPPP,KAPX) 
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40  WAY(KA)=WWAPEX(KNPPP,KAPX+1) 

CALL  MAPAZM(NXMAP,NYMAP,XM(KFF),WAX,WAY,WPAR) 
45  CONTINUE 
RETURN 
END 
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